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Abstract: Neutrino oscillations involving eV-scale neutrino mass states are investigated 
in the context of global neutrino oscillation data including short and long-baseline ac- 
celerator, reactor, and radioactive source experiments, as well as atmospheric and solar 
neutrinos. We consider sterile neutrino mass schemes involving one or two mass-squared 
differences at the eV^ scale denoted by 3+1, 3+2, and 1+3+1. We discuss the hints for 
eV-scale neutrinos from u e disappearance (reactor and Gallium anomalies) and v v e. 
appearance (LSND and MiniBooNE) searches, and we present constraints on sterile neu- 
trino mixing from 'i^'^ and neutral-current disappearance data. An explanation of all hints 
in terms of oscillations suffers from severe tension between appearance and disappearance 
data. The best compatibility is obtained in the 1+3+1 scheme with a p-value of 0.2% and 
exceedingly worse compatibilities in the 3+1 and 3+2 schemes. 
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1 Introduction 

Huge progress has been made in the study of neutrino oscillations [1-4], and with the recent 
determination of the last unknown mixing angle ^13 [5-10] a clear first-order picture of the 
three-flavor lepton mixing matrix has emerged, see e.g. [11]. Besides those achievements 
there are some anomalies which cannot be explained within the three-flavor framework and 
which might point towards the existence of additional neutrino flavors (so-called sterile 
neutrinos) with masses at the eV scale: 
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• The LSND experiment [12] reports evidence for i?^ — ?■ transitions with E/L ^ 
1 eV^, where E and L are the neutrino energy and the distance between source and 
detector, respectively. 

• This effect is also searched for by the MiniBooNE experiment [13-17], which reports 
a yet unexplained event excess in the low-energy region of the electron neutrino and 
anti-neutrino event spectra. No significant excess is found at higher neutrino energies. 
Interpreting the data in terms of oscillations, parameter values consistent with the 
ones from LSND are obtained. 

• Radioactive source experiments at the Gallium solar neutrino experiments SAGE and 
GALLEX have obtained an event rate which is somewhat lower than expected. This 
effect can be explained by the hypothesis of fe disappearance due to oscillations with 
Am^ > 1 eV^ [18, 19] ("Gallium anomaly"). 

• A recent re-evaluation of the neutrino flux emitted by nuclear reactors [20, 21] has led 
to somewhat increased fluxes compared to previous calculations [22-25]. Based on 
the new flux calculation, the results of previous short-baseline (L < 100 m) reactor 
experiments are in tension with the prediction, a result which can be explained by 
assuming Vf, disappearance due to oscillations with Am^ ~ 1 eV^ [26] ("reactor 
anomaly" ) . 

Sterile neutrino oscillation schemes have been considered for a long time, see e.g. [27-30] 
for early references on four- neutrino scenarios. Effects of two sterile neutrinos at the eV 
scale have been considered first in [31, 32], oscillations with three sterile neutrinos have 
been investigated in [33, 34]. 

Thus, while the phenomenology of sterile neutrino models is well known, it has also 
been known for a long time that the LSND and MiniBooNE v e appearance signals are 
in tension with bounds from disappearance experiments [35-37], challenging an interpre- 
tation in terms of sterile neutrino oscillations. This problem remains severe, and in the 
following we will give a detailed discussion of the status of the u ^ ^ u e appearance hints 
from LSND and MiniBooNE in the light of recent global data. The situation is better for 
the hints for e disappearance from the reactor and Gallium anomalies, which are not in 
direct conflict with any other data. This somewhat ambiguous situation asks for an exper- 
imental answer, and indeed several projects are under preparation or under investigation, 
ranging from experiments with radioactive sources, short-baseline reactor experiments, to 
new accelerator facilities. A recent review on light sterile neutrinos including an overview 
on possible experimental tests can be found in [38]. 

In this paper we provide an extensive analysis of the present situation of sterile neu- 
trino scenarios. We discuss the possibility to explain the tentative positive signals from 
LSND and MiniBooNE, as well as the reactor and Gallium anomalies in terms of ster- 
ile neutrino oscillations in view of the global data. New ingredients with respect to our 
previous analysis [39] are the following. 

• We use latest data from the MiniBooNE 'u^ — t- Vg appearance searches [15-17]. Our 
MiniBooNE appearance analysis is now based on Monte Carlo events provided by the 



- 2 - 



collaboration taking into account realistic event reconstruction, correlation matrices, 
as well as oscillations of various background components in a consistent way. 

• We include the constraints on the appearance probability from E776 [40] and ICARUS [41]. 

• We include the Gallium anomaly in our fit. 

• We take into account constraints from solar neutrinos, the KamLAND reactor exper- 
iment, and LSND and KARMEN measurements of the reaction v^, + -"^^C —7-6" + ^^N. 

• The treatment of the reactor anomaly is improved and updated by taking into ac- 
count small changes in the predicted anti-neutrino fluxes as well as an improved 
consideration of systematic errors and their correlations. 

• We take into account charged-current (CC) and neutral-current (NC) data from the 
MINOS long-baseline experiment [42, 43]. 

• We include data on disappearance from MiniBooNE [44] as well as Vy, disappear- 
ance from a joint MiniBooNE/SciBooNE analysis [45]. 

• In our analysis of atmospheric neutrino data, we improve our formalism to fully take 
into account the mixing of Vf, with other active or sterile neutrino states. 

All the data used in this work are summarized in Tab. 1. For other recent sterile neutrino 
global fits see [34, 46, 47]. We are restricting our analysis to neutrino oscillation data; 
implications for kinematic neutrino mass measurements and neutrino-less double beta- 
decay data have been discussed recently in [48-50]. 

Sterile neutrinos at the eV scale also have implications for cosmology. If thermalized 
in the early Universe they contribute to the number of relativistic degrees of freedom (ef- 
fective number of neutrino species A^^eff)- A review with many references can be found 
in [38]. Indeed there might be some hints from cosmology for additional relativistic de- 
grees of freedoms (A^efF bigger than 3), coming mainly from CMB data, e.g. [47, 51-55]. 
Recently precise CMB data from the PLANCK satellite have been released [56]. Depend- 
ing on which additional cosmological data are used, Aeff values ranging from 3.30^q 5^ to 
3.62]']Q4g (uncertainties at 95% CL) are obtained [56]. Constraints from Big Bang Nucle- 
osynthesis on Acff have been considered recently in [57]. Apart from their contribution 
to A'eff, thermalized eV-scale neutrinos would also give a large contribution to the sum of 
neutrino masses, which is constrained to be below around 0.5 eV. The exact constraint 
depends on which cosmological data sets are used, but the most important observables are 
those related to galaxy clustering [51-54]. In the standard ACDM cosmology framework 
the bound on the sum of neutrino masses is in tension with the masses required to explain 
the aforementioned terrestrial hints [54]. The question to what extent such sterile neu- 
trino scenarios are disfavored by cosmology and how far one would need to deviate from 
the ACDM model in order to accommodate them remains under discussion [47, 58, 59]. 
We will not include any information from cosmology explicitly in our numerical analysis. 
However, we will keep in mind that neutrino masses in excess of few eV may become more 
and more difficult to reconcile with cosmological observations. 
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Experiment 


dof 


channel 


comments 


Short-baseline reactors 


76 


- 




SBL 


Long-baseline reactors 


39 


- 


— ^ 


LBL 


KamLAND 


17 


TJ 






Gallium 


4 


- 




SBL 


Solar neutrinos 


261 




— ^ 


Ue + NC data 


LSND/KARMEN ^^C 


32 


Ve - 


— )■ I^e 


SBL 


CDHS 


15 


V,, - 


— ^ l/^ 


SBL 


MiniBooNE v 


15 


V,, - 


Up 


SBL 


MiniBooNE v 


42 




Uu 
p 


SBL 


MINOS CC 


20 


Vp. - 


Up 


LBL 


MINOS NC 


20 


Vf, - 


Us 


LBL 


Atmospheric neutrinos 


80 


(— ) 

y ^, - 


{— ] 
-5> p 


+ NC matter effect 


LSND 


11 


v^, - 


Ue 


SBL 


KARMEN 


9 


Vp. - 


Ue 


SBL 


NOMAD 


1 


Vp, - 


i^e 


SBL 


MiniBooNE v 


11 


Vp, ^ 


~5> J^e 


SBL 


MiniBooNE v 


11 


Vp - 




SBL 


E776 


24 


(— ) 
Vp - 


(— ) 


SBL 


ICARUS 


1 


Vp - 




LBL 


total 689 



Table 1. Summary of the data used in this work divided into 'z^'e, disappearance, and ap- 
pearance data. The column "dof" gives the number of data points used in our analysis minus the 
number of free nuisance parameters for each experiment. 

The outline of the paper is as follows. In Sec. 2 we introduce the formalism of sterile 
neutrino oscillations and fix the parametrization of the mixing matrix. We then consider 
e disappearance data in Sec. 3, discussing the reactor and Gallium anomalies. Constraints 
from ^ disappearance as well as neutral-current data are discussed in Sec. 4, and global 

— )• Ve appearance data including the LSND and MiniBooNE signals in Sec. 5. The 
global fit of all these data combined is presented in Sec. 6 for scenarios with one or two 
sterile neutrinos. We summarize our results and conclude in Sec. 7. Supplementary material 
is provided in the appendices including a discussion of complex phases in sterile neutrino 
oscillations, oscillation probabilities for solar and atmospheric neutrinos, as well as technical 
details of our experiment simulations. 

2 Oscillation parameters in the presence of sterile neutrinos 

In this work we consider the presence of s = 1 or 2 additional neutrino states with masses 
in the < few eV range. When moving from 1 to 2 sterile neutrinos the qualitative new 
feature is the possibility of CP violation already at short-baseline [33, 60].^ The neutrino 
mass eigenstates ui, . . . , 2/3+^ are labeled such that z^i , 1/2 , f 3 contribute mostly to the 
active flavor eigenstates and provide the mass squared differences required for "standard" 
three-flavor oscillations, Am^i ~ 7.5 x 10~^ eV^ and lAmg^^l « 2.4 x 10~^ eV^, where 

^Adding more than two sterile neutrinos does not lead to any qualitatively new physical effects and as 
shown in [33] the fit does not improve significantly. Therefore, we restrict the present analysis to s < 2 
sterile neutrinos. 
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Amjj = mf — rrij. The mass states U4, are mostly sterile and provide mass-squared 
differences in the range 0.1 eV^ < |Am|i|, lAmiil < 10 eV^. In the case of only one 
sterile neutrino, denoted by "3+1" in the following, we always assume Am|]^ > 0, but the 
oscillation phenomenology for Am|]^ < would be the same. For two sterile neutrinos, we 
distinguish between a mass spectrum where Aml;^ and Amg;^ are both positive ("3+2") 
and where one of them is negative ("1+3+1"). The phenomenology is slightly different 
in the two cases [61]. We assume that the s linear combinations of mass states which are 
orthogonal to the three flavor states participating in weak interactions are true singlets and 
have no interaction with Standard Model particles. Oscillation physics is then described by 
a rectangular mixing matrix \Jai with a = e, /u, r and i = I, 3 + s, and U*^Uf3i = Sap-"^ 

We give here expressions for the oscillation probabilities in vacuum, focusing on the 
3+2 case. It is trivial to recover the 3+1 formulas from them by simply dropping all 
terms involving the index "5". Formulas for the 1+3+1 scenario are obtained by taking 
either Am|^ or Am|]^ negative. Oscillation probabilities relevant for solar and atmospheric 
neutrinos are given in appendices C and D, respectively. 

First we consider the so-called "short-baseline" (SBL) limit, where the relevant range 
of neutrino energies and baselines is such that effects of Am|j^ and Am|]^ can be neglected. 
Then, oscillation probabilities depend only on Am?^ and Uai with i > 4. We obtain for 
the appearance probability 

P^^h'up^ = 4 |C/„4| sin^ ,^41 + 4 |C/,5| sin2 ^51 

+ S\UaAUi3iUa5Up^\ sin (/)4i sin 051 COs(</>54 - 7a/3) , (2.1) 

with the definitions 

Am?L 

(t^ij = , 7a/3 = arg (/a/354) , la/Bij = U*iUf3iUajU^j ■ (2.2) 

Eq. (2.1) holds for neutrinos; for anti-neutrinos one has to replace ^ai3 — 7a/3- Since 
Eq. (2.1) is invariant under the transformation 4 -f-)- 5 and 'jap —la/St we can restrict the 
parameter range to Atti^^ > 0, or equivalently Am^]^ > Am|^, without loss of generality. 
Note also that the probability Eq. (2.1) depends only on the combinations |f7a4f^/34| and 
\Ua5U[^^\. The only SBL appearance experiments we are considering are in the — )• 
V'e channel. Therefore, the total number of independent parameters is 5 if only SBL 
appearance experiments are considered. 

The 3+2 survival probability, on the other hand, is given in the SBL approximation 

by 

= 1 - 4 I 1 - 5^ iC/ailM ^ If^ml' sin^ cPa - 4 |?7„4| Vasi' sin' '^54 . (2.3) 

y i=4,5 / i=4,5 

In this work we include also experiments for which the SBL approximation cannot be 
adopted, in particular MINOS and ICARUS. For these experiments 03i is of order one. In 

^In this work we consider so-called phenomenological sterile neutrino models, where the 3 + s neutrino 
mass eigenvalues and the mixing parameters Ud are considered to be completely independent. In particular 
we do not assume a seesaw scenario, where the Dirac and Majorana mass matrices of the sterile neutrinos are 
the only source of neutrino mass and mixing. For such "minimal" sterile neutrino models see e.g. [62-64]. 
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the following we give the relevant oscillation probabilities in the limit of (j)4^i , (p^i , (p^^ — >■ oo 
and (j)2i — )• 0. We call this the long-baseline (LBL) approximation. In this case we obtain 
for the neutrino appearance probability (a / /3) 

5 

= 4|C/,3|'|f//33|' Sin2 (/,3i +Y^\U^,\^\UfSi\^ + 2Re(/„/345) 

4=4 

+ 4Re(lQ,/343 + /a/353) sin^ '^31 - 2Im(/Q,/343 + /a/353) sin(2(/)3i) . (2.4) 

The corresponding expression for anti-neutrinos is obtained by the replacement la/Sij 
I^pij- The survival probability in the LBL limit can be written as 

P};^^^^:'=[^-jZ\Uo.^?] +Elf^-l' + 2fl-Elf^-l')|f^"3|'cos(2</,3i). (2.5) 

\ i=3 / i=3 \ 1=3 / 

Note that in the numerical analysis of MINOS data neither the SBL nor the LBL approx- 
imations can be used because (/)3i, </)4i and can all become of order one either at the 
far detector or at the near detector [65]. Moreover, matter effects cannot be neglected in 
MINOS. All of these effects are properly included in our numerical analysis of the MINOS 
experiment. 

Sometimes it is convenient to complete the 3 x (3 -|- s) rectangular mixing matrix by 
s rows to an n X n unitary matrix, with n = 3 -|- s. For n = 5 we use the following 
parametrization for U: 

U = V^,>OMV2^V2A02zOi^OiiVnVi2 (2.6) 

where Oij represents a real rotation matrix by an angle 6ij in the ij plane, and Vij represents 
a complex rotation by an angle Oij and a phase (pij. The particular ordering of the rotation 
matrices is an arbitrary convention which, however, turns out to be convenient for practical 
reasons.^ We have dropped the unobservable rotation matrix V45 which just mixes sterile 
states. There is also some freedom regarding which phases are removed by field redefinitions 
and which ones are kept as physical phases. In appendix A we give a specific recipe for 
how to remove unphysical phases in a consistent way. Throughout this work we consider 
only phases which are phenomenologically relevant in neutrino oscillations. Under certain 
approximations, more phases may become unphysical. For instance, if an angle which 
corresponds to a rotation which can be chosen to be complex is zero the corresponding 
phase disappears. In practical situations often one or more of the mass-squared differences 
can be considered to be zero, which again implies that some of the angles and phases will 
become unphysical. In Tab. 2 we show the angle and phase counting for the SBL and LBL 
approximations for the 3-1-2 and 3-1-1 cases. 

In the notation of Eqs. (2.1), (2.3), (2.4), (2.5), it is explicit that only appearance 
experiments depend on complex phases in a parametrization independent way. However, in 
a particular parametrization such as Eq. (2.6), also the moduli \ Uai\ may depend on cosines 
of the phase parameters (pij , leading to some sensitivity of disappearance experiments to the 

^Note that the ordering chosen in Eq. (2.6) is equivalent to U = V35V25O15O34V24O14O23 V'i3Vi2, where 
the standard three-flavor convention appears to the right (apart from an additional complex phase), and 
mixing involving the mass states V4, and appear successively to the left of it. 
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A/P 


LBL approx. 


(A/P) 


SBL approx. 


(A/P) 


3+2 9/5 


V35V34V25O24O23O15O14V13 


(8/4) 


^35034^25024015014 


(6/2) 


3+1 6/3 


V34O24O23O14V13 


(5/2) 


O34O24O14 


(3/0) 



Table 2. Mixing angle and Phase counting for s = 2 (3+2) and s = 1 (3+1) sterile neutrino 
schemes. The column "A/P" denotes the number of physical angles and phases, respectively. The 
column "LBL approx." ("SBL approx.") corresponds to the approximation Am^i — (Amji — > 0, 
Amg]^ — > 0). We give also specific examples for which angles can be chosen real, by denoting with 
Vij (Oij) a complex (real) rotation. 



ipij in a CP-even fashion. Our parametrization Eq. (2.6) guarantees that 'z^'e disappearance 
experiments are independent of ipij. 

3 z/e and z/g disappearance searches 

Disappearance experiments in the z^'e sector probe the moduli of the entries in the first 
row of the neutrino mixing matrix, \Uei\- In the short-baseline limit of the 3+1 scenario, 
the only relevant parameter is \Ue4\- For two sterile neutrinos, also |C/e5| is relevant. In this 
section we focus on 3+1 models, and comment only briefly on 3+2. For 3+1 oscillations 
in the SBL limit, the z^e survival probability takes an effective two flavor form 

where we have defined an effective Vg-disappearance mixing angle by 

sin2 20ge = 4|C/e4p(l-|C/e4p). (3.2) 

This definition is parametrization independent. Using the specific parametrization of 
Eq. (2.6) it turns out that 9ee = On- 

3.1 SBL reactor experiments 

The data from reactor experiments used in our analysis are summarized in Tab. 3. Our 
simulations make use of a dedicated reactor code based on previous publications, see 
e.g. [79, 80]. We have updated the code to include the latest data and improved the 
treatment of uncertainties, see appendix B for details. The code used here is very similar 
to the one from Ref. [11], extended to sterile neutrino oscillations. The reactor experi- 
ments listed in Tab. 3 can be divided into short-baseline (SBL) experiments with baselines 
< 100 m, long-baseline (LBL) experiments with 100 m < L < 2 km, and the very long- 
baseline experiment KamLAND with an average baseline of 180 km. SBL experiments are 
not sensitive to standard three- flavor oscillations, but can observe oscillatory behavior for 
Aml]^, Amg]^ ~ 1 eV^. On the other hand, for long-baseline experiments, oscillations due 
to (Am|^, ^13) are most relevant, and oscillations due to eV^-scale mass-squared differences 
are averaged out and lead only to a constant flux suppression. KamLAND is sensitive to 
oscillations driven by (Am^i and ^12), whereas all 6ik with /c > 3 lead only to a constant 
flux reduction. 
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17 bins 





Table 3. Reactor data used in our analysis. The experiments in the upper part of the table 
have baselines L < 100 m and are referred to as SBL reactor experiments. For these experiments 
we list the baseline, the ratio of the observed and predicted rates (based on the flux predictions 
from [20, 21]), the uncorrelated error, and the total experimental error (i.e., the square-root of the 
diagonal entry of the correlation matrix) . Uncertainties from the neutrino flux prediction are not 
included here, but are taken into account in our numerical analysis. For details on the correlations 
and flux errors see appendix B. In the lower part of the table, we list experiments with baselines of 
order 1 km (LBL reactors), and the KamLAND experiment with an average baseline of 180 km. For 
DayaBay, RENO, and KamLAND, we do not give a number for the baseline here because several 
baselines are involved in each of these experiments. The number of SBL data points is 19 or 76 and 
the total number of reactor data points is 75 or 132, depending on whether a total rates analysis 
(3 data points) or a spectral analysis (25+25+10 bins) is used for the Bugey3 experiment. 

For the SBL reactor experiments we show in Tab. 3 also the ratio of the observed and 
predicted rate, where the latter is based on the flux calculations of [21] for neutrinos from 
239p^^ ^^^Pu fission and [20] for 238^ fission. The ratios are taken from [38] (which 
provides and update of [26]) and are based on the Particle Data Group's 2011 value for 
the neutron lifetime, = 881.5 s [81].^ We observe that most of the ratios are smaller 
than one. In order to asses the significance of this deviation, a careful error analysis is 
necessary. In the last column of Tab. 3, we give the uncorrelated errors on the rates. They 

*This number differs from their 2012 value 880.1 s by less than 0.2% [82]. The neutron lifetime enters 
the calculation of the detection cross section and therefore has a direct impact on the expected rate. The 
quoted uncertainties of < 0.2% are small compared to the uncertainties on the predicted flux, see [83] for 
a discussion of the neutron lifetime determination. 
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sm^29u Am|i[eV2] xLn/dof (GOF) A^Lcc/dof (CL) 



SBL rates only 


0.13 


0.44 


11.5/17 (83%) 


11.4/2 (99.7%) 


SBL incl. BugeyS spectr. 


0.10 


1.75 


58.3/74 (91%) 


9.0/2 (98.9%) 


SBL + Gallium 


0.11 


1.80 


64.0/78 (87%) 


14.0/2 (99.9%) 


SBL + LBL 


0.09 


1.78 


93.0/113 (92%) 


9.2/2 (99.0%) 


global i^e disapp. 


0.09 


1.78 


403.3/427 (79%) 


12.6/2 (99.8%) 



Table 4. Best fit oscillation parameters and xLin values as well as Axno-osc = Xno-osc ~ Xmin within 
a 3+1 framework. Except in the row labeled "SBL rates only", we always include spectral data 
from BugeyS. The row "global Ve disapp." includes the data from reactor experiments (see Tab. 3) 
as well as Gallium data, solar neutrinos and the LSND/KARMEN I'e disappearance data from 
j/g-^^C scattering. The CL for the exclusion of the no oscillation hypothesis is calculated assuming 
2 degrees of freedom (|J7e4| and Am|^). 

include statistical as well as uncorrelated experimental errors. In addition to these, there 
are also correlated experimental errors between various data points which are described in 
detail in appendix B. Furthermore, we take into account the uncertainty on the neutrino 
flux prediction following the prescription given in [21], see also appendix B for details. 

Fitting the SBL data to the predicted rates we obtain xVdof = 23.0/19 which corre- 
sponds to a p-value of 2.4%. When expressed in terms of an energy-independent normal- 
ization factor /, the best fit is obtained at 

/ = 0.935 ±0.024, Xmin/dof= 15.7/18 (p = 61%) , Ax/=i = 7.25 {2.7a). (3.3) 

Here Axj^^ denotes the improvement in compared to a fit with / = 1. Clearly the 
p- value increases drastically when / is allowed to float, leading to a preference for / 7^ 1 at 
the 2.7cj confldence level. This is our result for the significance of the reactor anomaly. Let 
us mention that (obviously) this result depends on the assumed systematic errors. While 
we have no particular reason to doubt any of the quoted errors, we have checked that 
when an adhoc additional normalization uncertainty of 2% (3%) is added, the significance 
is reduced to 2.1a {1.7a). This shows that the reactor anomaly relies on the control of 
systematic errors at the percent level. 

The flux reduction suggested by the reactor anomaly can be explained by sterile neu- 
trino oscillations. In Tab. 4 we give the best flt points and values obtained by fitting 
SBL reactor data in a 3+1 framework. The allowed regions in Am|]^ and sin 20\^ are 
shown in Fig. 1 (left) for a rate-only analysis as well as a fit including also Bugey3 spec- 
tral data. Both analyses give consistent results, with the main difference being that the 
spectral data disfavors certain values of Aml^^ around 0.6 — 0.7 eV^ and 1.3 eV^. The 
right panel of Fig. 1 shows the predicted rate suppression as a function of the baseline 
compared to the data. We show the prediction for the two best fit points from the left 
panel as well as one point located in the island around Am|^ ~ 0.9 eV^, which will be 
important in the combined fit with SBL appearance data. We observe that for the rate- 
only best fit point with Amli = 0.44 eV^ the prediction follows the tendency suggested 
by the ILL, Bugey4, and SRP (24 km) data points. This feature is no longer present for 
Amli > 1 eV^, somewhat preferred by Bugey3 spectral data, where oscillations happen at 
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Figure 1. Left: Allowed regions of oscillation parameters from SBL reactor data in the 3+1 
scheme for a rates only analysis (contours) as well as a fit including Bugey3 spectral data (colored 
regions). Right: Event rates in SBL reactor experiments compared to the predictions for three 
representative sets of oscillation parameters. The thick (thin) error bars correspond to uncorrelated 
(total) experimental errors. The neutrino flux uncertainty is not included in the error bars. The 
Rovno and SRP data points at 18 m have been shifted for better visibility. 



even shorter baselines. However, from the GOF values given in Tab. 4 we conclude that 
also those solutions provide a good fit to the data. 

3.2 The Gallium anomaly 

The response of Gallium solar neutrino experiments has been tested by deploying radioac- 
tive ^^Cr or 37 Ar sources in the GALLEX [84, 85] and SAGE [86, 87] detectors. Results 
are reported as ratios of observed to expected rates, where the latter are traditionally 
computed using the best fit cross section from Bahcall [88], see e.g. [19]. The values for 
the cross sections weighted over the 4 (2) neutrino energy lines from Cr (Ar) from [88] 
are cJB(Cr) = 58.1 x 10~^^cm^, (TB(Ar) = 70.0 x 10"'^^ cm^. While the cross section for 
7^Ga — )• ''^Ge into the ground state of ^^Ge is well known from the inverse reaction there 
are large uncertainties when the process proceeds via excited states of "^^Ge at 175 and 
500 keV. Following [88], the total cross section can be written as 

fv\ fv\f'\ I BGT175 BGT5oo\ , .s 

a{X) = a,.,.{X) ^1 + ax^^ + bx^^J (3.4) 

with X = Cr, Ar. The coefficients ocr = 0.669, bcr = 0.220, OAr = 0.695, &Ar = 0.263 are 
phase space factors. The ground state cross sections are precisely known [88]: (Tg.s.(Cr) = 
55.2 X 10-^6 cm2, crg.s.(Ar) = 66.2 x lO'^^cm^. BGT denote the Gamov-Teller strength of 
the transitions, which have been determined recently by dedicated measurements [89] as 

1^^^ = 0.0399 ± 0.0305 , = 0.207 ± 0.016 . (3.5) 

BGTg.s. BGTg.s. 

In our analysis we use these values together with Eq. (3.4) for the cross section. 
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Aml^ 


^14 


^15 


xLn (GOF) 


Axi+i (CL) 


AxLosc (CL) 


SBLR 


0.46 


0.87 


0.12 


0.13 


53.0/(76-4) (95%) 


5.3 (93%) 


14.3 (99.3%) 


SBLR+gal 


0.46 


0.87 


0.12 


0.14 


60.2/(80-4) (90%) 


3.8 (85%) 


17.8 (99.9%) 



Table 5. Best fit point of SBL reactor data and SBL reactor + Gallium data in a 3+2 oscillation 
scheme. We give the mass-squared differences in eV^ and the mixing angles in radians. The relation 
to the mixing matrix elements is \Ue4\ — cos ^15 sin 6*14 and \Ue5\ = sin6'i5. The Ax^ relative to 
3+1 oscillations is evaluated for 2 dof, corresponding to the two additional parameters, while for 
the Ax^ relative to no-oscillations we use 4 dof. 



This means that the ratios of observed to expected rates based on the Bahcall predic- 
tion have to be rescaled by a factor 0.982 (0.977) for the Cr (Ar) experiments, so that we 
obtain for them the following updated nmnbers for our fits: 

i?i(Cr) = 0.94 ± 0.11 [85] ^^^^ ( i?3(Cr) = 0.93 ± 0.12 [86] 
i?2(Cr) = 0.80 ± 0.10 [85] ' ' \ i?4(Ar) = 0.77 ± 0.08 [87] " 

Here, we have symmetrized the errors, and we have included only experimental errors, but 
not the uncertainty on the cross section (see below). 

We build a out of the four data points from GALLEX and SAGE and introduce two 
pulls corresponding to the systematic uncertainty of the two transitions to excited state 
according to Eq. (3.5). The determination of BGT175 is relatively poor, with zero being 
allowed at 2a. In order to avoid unphysical negative contributions from the 175 keV state, 
we restrict the domain of the corresponding pull parameter accordingly. Fitting the four 
data points with a constant neutrino flux normalization factor r we find 

xLn = 2.26/3 dof, r^i„ = 0.84t0:°54 , Ax,Li = 8.72 (2.95a) (3.7) 

Because of the different cross sections used, these results differ from the ones obtained 
in [19], where the best fit point is at r = 0.76, while the significance is comparable, around 
3(7. An updated analysis including also a discussion of the implications of the measurement 
in [89] can be found in [90] . 

The event deficit in radioactive source experiments can be explained by assuming 
mixing with an eV-scale state, such that Ue disappearance happens within the detector 
volume [18]. We fit the Gallium data in the 3+1 framework by averaging the oscillation 
probability over the detector volume using the geometries given in [18]. The resulting 
allowed region at 95% confidence level is shown in orange in Fig. 2. Consistent with 
the above discussion we find mixing angles somewhat smaller than those obtained by the 
authors of [19]. The best fit point from combined Gallium+SBL reactor data is given 
in Tab. 4, and the no-oscillation hypothesis is disfavored at 99.9% CL (2 dof) or 3.3a 
compared to the 3+1 best fit point. 

Let us consider now the Gallium and SBL reactor data in the framework of two sterile 
neutrinos, in particular in the 3+2 scheme. SBL and z^e disappearance data depend 
on 4 parameters in this case, Am-I,, Am-5,, and the two mixing angles ^14 and ^15 (or, 
equivalently, the moduli of the two matrix elements Uei and C/es)- We report the best 
fit points from SBL reactor data and from SBL reactor data combined with the Gallium 
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Figure 2. Allowed regions at 95% CL (2 dof) for 3+1 oscillations. We show SBL reactor data 
(blue shaded), Gallium radioactive source data (orange shaded), disappearance constraints 
from Vf.-^'^C scattering data from LSND and KARMEN (dark red dotted), long-baseline reac- 
tor data from CHOOZ, Palo Verde, DoubleChooz, Daya Bay and RENO (blue short-dashed) and 
solar+KamLAND data (black long-dashed). The red shaded region is the combined region from all 
these Ve and v^. disappearance data sets. 

source data in Tab. 5. For these two cases we find an improvement of 5.3 and 3.8 units in 
X^, respectively, when going from the 3-1-1 scenario to the 3-1-2 case. Considering that the 
2>+2 model has two additional parameters compared to 3-1-1, we conclude that there is no 
improvement of the fit beyond the one expected by increasing the number of parameters, 
and that SBL v ^ data sets show no significant preference for 3-1-2 over 3+1. This is 
also visible from the fact that the confidence level at which the no oscillation hypothesis is 
excluded does not increase for 3+2 compared to 3+1, see the last columns of Tabs. 4 and 5. 
There the Ax^ is translated into a confidence level by taking into account the number of 
parameters relevant in each model, i.e., 2 for 3+1 and 4 for 3+2. 

3.3 Global data on Vf. and Vf. disappearance 

Let us now consider the global picture regarding disappearance. In addition to the 
short-baseline reactor and Gallium data discussed above, we now add data from the fol- 
lowing experiments: 

• The remaining reactor experiments at a long baseline ( "LBL reactors" ) and the very 
long-baseline reactor experiment KamLAND, see table 3. 

• Global data on solar neutrinos, see appendix C for details. 

• LSND and KARMEN measurements of the reaction i/e + ^^C e~ + ^^N [91, 92]. 
The experiments have found agreement with the expected cross section [93], hence 
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their measurements constrain the disappearance of with eV-scale mass-squared 
differences [94, 95]. Details on our analysis of the ^^C scattering data are given in 
appendix E.l. 

So far the LBL experiments DayaBay and RENO have released only data on the rela- 
tive comparison of near (L ~ 400 m) and far (L ~ 1.5 km) detectors, but no information on 
the absolute flux determination is available. Therefore, their published data are essentially 
insensitive to oscillations with eV-scale neutrinos and they contribute only indirectly via 
constraining ^13. In our analysis we include a free, independent flux normalization factor 
for each of those two experiments. Chooz and DoubleChooz both lack a near detector. 
Therefore, in the official analyses performed by the respective collaborations the Bugey4 
measurement is used to normalize the flux. This makes the official Chooz and DoubleChooz 
results on 6*13 also largely independent of the presence of sterile neutrinos. However, the 
absolute rate of Bugey4 in terms of the flux predictions is published (see Tab. 3) and we 
can use this number to obtain an absolute flux prediction for Chooz and DoubleChooz. 
Therefore, in our analysis Chooz and DoubleChooz (as well as Palo Verde) by themselves 
also show some sensitivity to sterile neutrino oscillations. In a combined analysis of Chooz 
and DoubleChooz with SBLR data the official analyses are recovered approximately. Pre- 
vious considerations of LBL reactor experiments in the context of sterile neutrinos can be 
found in Refs. [96-99]. 

We show in Tab. 4 a combined analysis of the SBL and LBL reactor experiments 
(row denoted by "SBL+LBL"), where we minimize with respect to ^13. We find that the 
significance of the reactor anomaly is not affected by the inclusion of LBL experiments and 
finite ^13. The Axno-osc even slightly increases from 9.0 to 9.2 when adding LBL data to the 
SBL data ("no-osc" refers here to ^14 = 0). Hence, we do not agree with the conclusions 
of Ref. [100], which finds that the significance of the reactor anomaly is reduced to \Aa 
when LBL data and a finite value of ^13 is taken into account. 

Solar neutrinos are also sensitive to sterile neutrino mixing (see e.g. [101-103]). The 
main effect of the presence of Vf. mixing with eV states is an over-all flux reduction. While 
this effect is largely degenerate with ^13, a non-trivial bound is obtained in the combination 
with DayaBay, RENO and KamLAND. KamLAND is sensitive to oscillations driven by 
Am^i and ^12, whereas sterile neutrinos affect the overall normalization, degenerate with 
^13. The matter effect in the sun as well as SNO NC data provide additional signatures 
of sterile neutrinos, beyond an overall normalization. As we will show in Sec. 4 solar data 
depend also on the mixing angles ^24 and ^34, controlling the fraction of Vf. — )• transitions, 
see e.g. [101]. As discussed in appendix C, in the limit Am?^ = 00 for i > 3, solar data 
depends on 6 real mixing parameters, 1 complex phase and Aml^. Hence, in a 3+1 scheme 
all six mixing angles are necessary to describe solar data in full generality. However, once 
other constraints on mixing angles are taken into account the effect of ^24, ^34, and the 
complex phase are tiny and numerically have a negligible impact on our results. Therefore 
we set ^24 = ^34 = for the solar neutrino analysis in this section. In this limit solar data 
becomes also independent of the complex phase. 

The results of our fit to global v f. disappearance data are shown in Fig. 2 and the 
best fit point is given in Tab. 4. For this analysis the mass-squared differences Am2i and 
Arrig]^ have been fixed, whereas we marginalize over the mixing angles 612 and ^13. We see 
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|t/e4|' |t/e4|' 

Figure 3. Constraints on and disappearance in a 3 + 1 model at two different fixed values of 
Arn|]^. Regions are shown at 95% CL (2 dof) with respect to the minimum at the fixed /S.m\^. 
We show constraints from the radio chemical Gallium experiments using radioactive sources (orange 
band), from short-baseline reactor experiments (blue band), from the KARMEN and LSND mea- 
surements of the z^g-^^C cross section (dark red dotted line), from long-baseline reactor experiments 
(blue dashed line), from the combined solar -|-KamL AND data (black dashed line), and from the 
the combination of all aforementioned experiments (red shaded region). 



from Fig. 2 that the parameter region favored by short-baseline reactor and Gallium data 
is well consistent with constraints from long-baseline reactors, KARMEN's and LSND's Ve 
rate, and with solar and KamLAND data. 

Recently, data from the Mainz [104] and Troitsk [105] tritium beta-decay experiments 
have been re-analyzed to set limits on the mixing of with new > eV neutrino mass 
states. Taking the results of [105] at face value, the Troitsk limit would cut-off the high- 
mass region m Fig. 2 at around 100 eV^ [106] (above the plot -range shown in the figure). 
The bounds obtained in [104] are somewhat weaker. The differences between the limits 
obtained in [104] and [105] depend on assumptions concerning systematic uncertainties 
and therefore we prefer not to explicitly include them in our fit. The sensitivity of future 
tritium decay data from the KATRIN experiment has been estimated in [107]. Implications 
of sterile neutrinos for neutrino-less double beta-decay have been discussed recently in [48- 
50]. 

Let us now address the question whether the presence of a sterile neutrino affects the 
determination of the mixing angle ^13 (see also [99, 100]). In Fig. 3 we show the combined 
determination of ^13 and On for two fixed values of Amlj^ . The left panel corresponds to a 
relatively large value of 10 eV^, whereas for the right panel we have chosen the value favored 
by the global v e disappearance best fit point, 1.78 eV^. The mass-squared differences Am|^ 
and Am^^ have been fixed, whereas we marginalize over the mixing angle 9i2- We observe a 
clear complementarity of the different data sets: SBL reactor and Gallium data determine 
|C/e4|, since oscillations are possible only via Am|]^, all other mass-squared differences are 
effectively zero for them. For LBL reactors Am^^ can be set to infinity, ^rn\i is finite, 
and Am^]^ is effectively zero; therefore they provide an unambiguous determination of ^13 
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by comparing near and far detector data. The upper bound on |C/e4| from LBL reactors is 
provided by Chooz, Palo Verde, DoubleChooz, since for those experiments also information 
on the absolute flux normalization can be used, as mentioned above. In contrast, for solar 
neutrinos and KamLAND, both Amli and Arrig]^ are effectively infinite, and ^13 and 614 
affect essentially the overall normalization and are largely degenerate, as visible the figure. 

In conclusion, the ^13 determination is rather stable with respect to the presence of 
sterile neutrinos. We note, however, that its interpretation becomes slightly more com- 
plicated. For instance, in the 3+1 scheme using the parametrization from Tab. 2, the 
relation between mixing matrix elements and mixing angles is | Ues \ = cos ^14 sin ^13 and 
\Ue4,\ = sin ^14. Hence, the one-to-one correspondence between {Uesl in the 

three-flavor case is spoiled. 

4 z/^, u^, and neutral-current disappearance searches 

In this section we discuss the constraints on the mixing of and with new eV-scale 
mass states. In the 3+1 scheme this is parametrized by \Ufj,i\ and {Ur^l, respectively. 
In terms of the mixing angles as defined in Eq. (2.6) we have |C//x4| = cos0i4sin024 and 
\Ur4\ = cos ^14 cos 024 sin ^34. In the present paper we include data sets from the following 
experiments to constrain and i^V mixing with eV states: 

• SBL disappearance data from CDHS [108]. Details of our simulation are given 
in [79]. 

• Super-Kamiokande. It has been pointed out in [109] that atmospheric neutrino data 
from Super-Kamiokande provide a bound on the mixing of with eV-scale mass 
states, i.e., on the mixing matrix elements |C//^4|, [f^/^sl- In addition, neutral-current 
matter effects provide a constraint on |?7t4| , l^^rsl • A discussion of the effect is given in 
the appendix of [33]. Details on our analysis and references are given in appendix D. 

• MiniBooNE [44, 45]. Apart from the *z^'e appearance search, MiniBooNE can also 
look for SBL disappearance. Details on our analysis are given in appendix E.4. 

• MINOS [42, 43]. The MINOS long-baseline experiment has published data on charged 
current (CC) disappearance as well as on the neutral current (NC) count rate. 
Both are based on a comparison of near and far detector measurements. In addition 
to providing the most precise determination of Am|^ (from CC data), those data can 
also be used to constrain sterile neutrino mixing, where CC (NC) data are mainly 
relevant for \U^4\, {U^^l {\Uri\, [C^tsD- See appendix E.5 for details. 

Additional constrains on mixing with eV-scale states (not used in this analysis) can be 
obtained from data from the Ice Cube neutrino telescope [110-115]. 

Limits on the row of the mixing matrix come from disappearance experiments. 
In a 3+1 scheme the SBL disappearance probability is given by 

pSBL,3+i ^ 1 _ _ ^.^2 Arr|L ^ ^ _ ^^^^ Arn^ ^ 
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Figure 4. Left: Constraints in the plane of |C/;i4p and Amli at 99% CL (2 dof) from CDHS, 
atmospheric neutrinos, MiniBooNE disappearance, MINOS CC and NC data, and the combination 
of them. We minimize with respect to \Ur4\ and the complex phase 9324 • In red we show the region 
preferred by LSND and MiniBooNE appearance data combined with reactor and Gallium data on 
V'e disappearance, where for fixed |J7^4p we minimize with respect to |J7e4p- Right: Constraints 
in the plane of \Uri\^ and Amj^ at 99% CL (2 dof) from MINOS CC + NC data (green) and the 
combined global and NC disappearance data (blue region, black curves). We minimize with 
respect to |C/^4| and we show the weakest ("best phase") and strongest ("worst phase") limits, 
depending on the choice of the complex phase (/524- In both panels we minimize with respect to 
Amf^, 023, and we fix sin^ 26*13 0.092 and 6*14 — (except for the evidence regions in the left 
panel). 

where we have defined an effective disappearance mixing angle by 

sin2 20^^ = 4|C/^4|'(1 - |C/mI') > (4-2) 

i.e., in our parametrization (2.6) the effective mixing angle 6*^^ depends on both 624 and 614. 
In contrast to the z^e disappearance searches discussed in the previous section, experiments 
probing ^ disappearance have not reported any hints for a positive signal. We show 
the limits from the data listed above in the left panel of Fig. 4. Note that the MINOS 
limit is based on the comparison of the data in near and far detectors. For Am|i ~ 
10 eV^ oscillation effects become relevant at the near detector, explaining the corresponding 
features in the MINOS bound around that value of Amli, whereas the features around 
Am|i ~ O.leV^ emerge from oscillation effects in the far detector. The roughly constant 
limit in the intermediate range 0.5 eV^ < Amli < 3eV^ corresponds to the limit Am|i ~ 
0(00) in the near (far) detector adopted in [42, 43]. In that range the MINOS limit on 
|f7^4| is comparable to the one from SuperK atmospheric data. For Am|^ > 1 eV^ the limit 
is dominated by CDHS and MiniBooNE disappearance data. 

In Fig. 4 (left) we show also the region preferred by the hints for eV-scale oscillations 
from LSND and MiniBooNE appearance data (see next section) combined with reactor 
and Gallium data on *j/'e disappearance. For fixed |[/^4p we minimize the corresponding 
with respect to |[/e4p to show the projection in the plane of |C//i4p and Aml^. The 



-16- 



tension between the hints in the — t- 'i^'g and Vg — ?■ channels compared to the limits 
from — )• 'j/'^ data is clearly visible in this plot. We will discuss this conflict in detail in 
section 6. 

Limits on the mixing of Ur with eV-scale states are obtained from data involving 
information from NC interactions, which allow to distinguish between z^^ — )• Vr and 
^ — )• s transitions.^ The relevant data samples are atmospheric and solar neutrinos (via 
the NC matter effect) and MINOS NC data. Furthermore, the parameter \Uri\ controls the 
relative weight of the oscillation modes — t- and — )• Vg at the "atmospheric" scale 
Aml]^: a large value of |f7r4| implies a large fraction of — )• oscillations at the Am|^ 
scale. The limit in the plane of |t/T4p and Am|^ is shown in the right panel of Fig. 4. 

As follows from Eq. (2.4) (see also appendix A), in the LBL approximation relevant for 
MINOS NC data a complex phase enters the oscillation probabilities, corresponding to the 
combination aig{U*^Ur4U^3U*^) . In our calculations we take the rotation matrix V24 to be 
complex and use the phase ip2i to parametrize this phase. In Fig. 4 we illustrate the impact 
of this phase by showing the strongest and weakest limit obtained when varying (^24- We 
observe that the limit from MINOS depends quite significantly on this phase. The different 
shapes of the "best phase" and "worst phase" regions emerge from the different properties 
of CC and NC data. For the weakest limit ("best phase") the fit uses the freedom of the 
term including the complex phase, which implies that a finite value of ^24 (or |C/^4|) is 
adopted, subject to the constraint from MINOS CC data. Therefore the same structure 
as in the left panel of Fig. 4 becomes visible also in limit on \Uri\- If we force the phase 
to take on a value not favored by the fit, a smaller is obtained for 6*24 close to zero, 
which implies that the phase actually becomes unphysical. In this case CC data are not 
important for the limit on |f/r4|) which then is dominated by NC data. Because of the much 
worse energy reconstruction for NC events compared to CC ones, the features induced by 
finite values of Am|^ in either the far or near detector become to a large extent washed 
out. 

The global limit on | Ur4 \ is actually dominated by atmospheric neutrino data and shows 
only a very weak dependence on the complex phase. In our atmospheric neutrino analysis 
the information on \Uta\ enters via the NC matter effect induced by the presence of sterile 
neutrinos. A large value of |?7r4| would imply a significant matter effect in Am|| driven 
'z/'^ disappearance, which is not consistent with the zenith angle distribution observed in 
SuperK. We find the limit 

|?7^4p<0.2 at 2cr(ldof) (4.3) 

from global data, largely independent of Am|^ as well as complex phases. 

Fig. 5 shows the constraints in the plane of |C//x4p and |C/T-4p for three fixed values of 
Aml;^. We observe the comparable bound on |f/^4p from MINOS (mainly CC data) and 
atmospheric, which however is superseded by CDHS, MiniBooNE for Am|]^ > 1 eV^ (left 
and middle panels). Those latter data however, do not provide any constraint on |f7,-4p, 
where the global bound is dominated by atmospheric neutrinos for all values of Am|^ of 
interest. We also observe that solar neutrinos provide a bound on |C/r4p of similar strength 

^The searches for Vr appearance at NOMAD [116] and CHORUS [117] at short baselines are sensitive 
only to specific parameter combinations like |(7^4t/T-4| or |C/e4'7T4| and therefore do not provide a constraint 
on |L''t4| by itself. 
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|f/;rf|' \U^\^ |!/p4|' 



Figure 5. Constraints in the plane of |J7/i4p and |f/r4p for three fixed values of Amli from MINOS 
CC + NC data (green) , atmospheric neutrinos (orange) , CDHS + MiniBooNE disappearance + 
LBL reactors (red), and the combination of those data (blue). The constraint from solar neutrinos is 
shown in magenta. Regions are shown at 90% and 99% CL (2 dof) with respect to the minimum 
at the fixed Aml-^ . We minimize with respect to complex phases and include effects of O13 and 614 
where relevant. The gray region is excluded by the unitarity requirement |C/^4p + |f/r4p < 1- Note 
the different scale on the axes. 



as MINOS data, thanks to the NC matter effect and SNO NC data. No relevant limit can 
be set on |?7^4p from solar neutrinos. 

5 z/^ — i- i/g and z/^ — )■ z/g appearance searches 

Now we move on to the discussion of appearance searches. In contrast to disappearance 
experiments which probe only one row of the mixing matrix, i.e., only the elements \Uai\ 
for fixed a, an appearance experiment in the channel Vq, — )• is sensitive to two rows via 
combinations like |C/Q,jC/gj| and potentially to some complex phases. In the SBL approxi- 
mation the 3+1 appearance probability in the phenomenologically most relevant channel 
ly ^ ^ v e takes the form 

pSBL,3+i ^ ^u,Ue,\' sin2 = W 2^ sm' , (5.1) 

where we have defined an effective mixing angle by 

Sm'^29f,e=MU^AUe4\^ ■ (5.2) 

In the parametrization from Eq. (2.6) we obtain sin20^e = sin 024 sin 2^14. The oscillation 
probability in the 3+2 scheme is given in Eq. (2.1). The 3+1 SBL appearance probability 
does not depend on complex phases, whereas in the 3+2 scheme CP violation via complex 
phases is possible at SBL [33, 60]. 

Our analyses of LSND [12], KARMEN [118], NOMAD [119] 'i^'^ appearance 
data are based on [33, 79, 120], where references and technical details can be found. Our 
analyses of E776 [40] and ICARUS [41], used for the first time in the present paper, are 
described in appendices E.2 and E.3, respectively.^ In the case of LSND, we use only the 



^Recently also the OPERA experiment presented results from a. ^ Ve appearance search [121]. The 
obtained limit is comparable to the one from ICARUS [41]. 
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decay-at-rest (DAR) data which are most sensitive to oscillations. Decay-in-flight (DIF) 
data on — )• Vf. are consistent with the signal seen in DAR data, however the significance 
of the oscillation signal for DIF is much less than for DAR. A combined DAR-DIF analysis 
in a two-neutrino framework would shift the allowed region to somewhat smaller values of 
the mixing angle. A detailed discussion of LSND DAR versus DIF in the context of 3+1 
neutrino oscillations can be found in [35]. 

In our analysis of the MiniBooNE and appearance search we use the latest data^ 
from [16], following closely the analysis instructions provided by the collaboration. Details 
are given in appendix E.4. Since their very first data release in 2007 [13], MiniBooNE 
observe an excess of events over expected background in the low energy (< 500 MeV) 
region of the event spectrum [122]. Since the spectral shape of the excess is difficult to 
explain in a two-flavor oscillation framework, historically the analysis window has been 
(somewhat artificially) divided into a low energy region containing the excess events and a 
high energy part with no excess.^ Preliminary results from anti-neutrinos showed also some 
indication for an event excess in the high energy part of the spectrum [14] which indicated 
the need for CP violation in order to reconcile neutrino and anti-neutrino data. However, 
for the most recent data [15, 16] the shapes of the neutrino and anti-neutrino spectra appear 
to be consistent with each other, showing excess events below around 500 MeV and data 
consistent with background in the high energy region, see Fig. 6. In our work we always 
analyse the full energy spectrum for both neutrinos and anti-neutrinos. Contrary to the 
analysis of the MiniBooNE collaboration we take into account oscillations of all background 
components in a consistent way, according to the particular oscillation framework to be 
tested, see appendix E.4 for details. 

In Fig. 7 we show a summary of the — )• v e data in the 3+1 scheme. We observe 
an allowed region from MiniBooNE anti-neutrino data that is driven by the event excess 
below around 800 MeV and has significant overlap with the parameter region preferred 
by LSND. At the 99% CL shown in the figure, MiniBooNE neutrino data give only an 
upper bound, although we find closed regions (again driven by the low-energy excess) at 
lower confidence levels. This is in qualitative agreement with the results obtained by the 
MiniBooNE collaboration, compare Fig. 4 of [16] or Fig. 3 of [17]. The different shape of 
our regions is due to the oscillations of the background components. Those can be relatively 
large in an appearance only fit, since for fixed sin^ 20^e we allow and |C/e4| to vary 

freely, subject to the constraint Eq. (5.2). We have checked that when we adopt the same 
assumptions as the MiniBooNE collaboration we recover their regions/bounds with good 
accuracy. 

The recent constraint on — )• appearance from ICARUS [41] at long-baseline 
leads to a bound on sin^ 20^e essentially independent of Am|^ in the range shown here. 
It excludes in particular the region of large mixing and low Am^]^ that is otherwise un- 
constrained by appearance experiments.^ An important background for the Amli driven 

^The recent updated analysis from MiniBooNE [17] is based on the same data as [16], corresponding to 
6.46 X 10^° protons on target in neutrino mode and 11.27 x 10^° protons on target in anti-neutrino mode. 

*The importance of energy reconstruction effects for tlie low energy excess has been pointed out in 
Refs. [123, 124], see also [17]. 

^Note that this region is also excluded by Ue and disappearance searches once Eq. (5.2) is used to 
relate sin^ 26^^ to the effective mixing angles probed by the disappearance experiments. 



-19- 




0.2 0.4 0.6 0.8 1. 1.2 1.4 1.5 3. 0.2 0.4 0.6 0.8 1. 1.2 1.4 1.5 3. 

(GeV) Ey (GeV) 



Figure 6. MiniBooNE neutrino (left) and anti-neutrino (right) data compared to the predicted 
spectra for the 3+1, 3+2, and 1+3+1 best fit points for the combined appearance data (the data 
set used in Fig. 7) and global data including disappearance. Shaded histograms correspond to the 
unoscillated backgrounds. The predicted spectra include the effect of background oscillations. The 
corresponding values (for combined neutrino and anti-neutrino data) are also given in the plot. 




Figure 7. Allowed regions and upper bounds at 99% CL (2 dof ) for — V'e appearance experi- 
ments in the 3+1 scheme. We show the regions from LSND and MiniBooNE anti-neutrino data and 
the bounds from MiniBooNE neutrinos, KARMEN, NOMAD, ICARUS, and E776. The latter is 
combined with LBL reactor data in order to constrain the oscillations of the V'e backgrounds; this 
leads to a non- vanishing bound on sin^ 20^2 from E776 at low Am|]^. The red region corresponds 
to the combination of those data, with the star indicating the best fit point. 

— )■ Vf, search in ICARUS are Ve appearance events due to Am|]^ and ^13. Furthermore, 
as discussed in section 2 and appendix A the long-baseline appearance probability in the 
3+1 scheme depends on one complex phase. In deriving the ICARUS bound shown in 



-20- 



Experiment 


X3+i/dof 


9 1 ^ r 

X3+2/dof 


XI+3+i/dof 


LSND 


11.0/11 


8.6/11 


7.5/11 




LiJ ,0 1 J- i 




9 1/11 


MiniB v 


10.7/11 


9.6/11 


12.7/11 
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0.0/1 


ICARUS 


2.0/1 


2.3/1 


1.5/1 


Combined 


87.9/(68 - 2) 


72.7/(68- 5) 


74.6/(68 - 5) 



Table 6. Individual contributions to the at the best fit point of the combined appearance data 
for 3+1, 3+2, and 1+3+1. The individual values do not add up to the number for the combined 
fit given in the last row because of correlations between MiniBooNE neutrino and anti-neutrino 
data. 

Fig. 7 we fix the parameters sin^ 2^13 = 0.092 and Am|^ = 2.4 x 10~^ eV^ but marginalize 
over the relevant complex phase. 

As visible in Fig. 7 there is a consistent overlap region for all v ^ ^ v e experiments 
and we can perform a combined analysis. The resulting region is shown in red in Fig. 7. The 
best fit point is at sin^ 26'^e = 0.013, Am|i = 0.42 eV^ with Xmin/dof = 87.9/(68 - 2) dof 
(GOF = 3.7%). The no-oscillation hypothesis is excluded with respect to the best fit point 
with Ax^ = 47.7. This large value is mostly driven by LSND. The relatively low GOF 
comes mainly from MiniBooNE neutrino data, as can be seen from Tab. 6, where we list the 
individual contribution of the experiments to the total appearance x^. This is also obvious 
from Fig. 6, showing that at the 3+1 appearance best fit point (black dotted histogram) 
the fit to the neutrino spectrum is not very good, predicting too much excess in the region 
0.6 — 1 GeV and only partially explaining the excess in the data below 0.4 GeV. 

Analysing the same data in the 3+2 scheme we find a best fit point at Am^j^ = 0.57 eV^, 
Am2^ = 1.24 eV^ with x^i^/dof = 72.7/(68 - 5) (GOF = 19%). The GOF improves 
considerably with respect to 3+1. We find 

X3+l,app ~ X3+2,app = 15.2 . (5.3) 

For 3 dof (corresponding to the 3 additional SBL appearance parameters in 3+2) this 
implies that appearance data favor 3+2 over 3+1 at the 99.8% CL. From Tab. 6 we see 
that basically all experiments have a reasonable x^/dof value (maybe with the exception 
of E776, which intrinsically has a somewhat high x^)- In particular MiniBooNE neutrino 
data improve by 8.7 units compared to 3+1. This is also visible in Fig. 6, with the red 
dotted curve (3+2 appearance best fit) showing a much better fit than the black dotted 
one (3+1 appearance best fit), with x^ = 24 for 22 dof for the joint MiniBooNE neutrino 
and anti-neutrino data. The appearance data fit in a 1+3+1 scheme is similar to the 3+2 
case, with a slightly better fit for LSND and MiniBooNE neutrino, and a slightly worse fit 
for MiniBooNE anti-neutrino data, compare Tab. 6. The predicted MiniBooNE spectra at 
the 1+3+1 appearance best fit are shown as blue dotted histograms in Fig. 6. We find for 
1+3+1 Xmin/dof= 74.6/(68 - 5) (GOF = 15%) and 

Xi+i,app - X?+3+i,app = 13.3 (99.6% CL, 3 dof). (5.4) 
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Figure 8. Results of the global fit in the 3+1 scenario, shown as exclusion limits and allowed regions 
for the effective mixing angle sin^ 20^e = 4|C/e4p|t^/i4p and the mass squared difference Am|]^. Left: 
Comparison of the parameter region preferred by appearance data (LSND, MiniBooNE appearance 
analysis, NOMAD, KARMEN, ICARUS, E776) to the exclusion hmit from disappearance data 
(atmospheric, solar, reactors. Gallium, CDHS, MINOS, MiniBooNE disappearance, KARMEN and 
LSND i/g-^^C scattering). Right: Regions preferred by experiments reporting a signal for sterile 
neutrinos (LSND, MiniBooNE, SBL reactors. Gallium) versus the constraints from all other data, 
shown separately for disappearance and appearance experiments, as well as their combination. 



6 Combined analysis of global data 

We now address the question whether the hints for sterile neutrino oscillations discussed 
above can be reconciled with each other as well as with all existing bounds within a com- 
mon sterile oscillation framework. In section 6.1 we discuss the 3+1 scenario, whereas in 
section 6.2 we investigate the 3+2 and 1+3+1 schemes. 

6.1 3+1 global analysis 

In the 3+1 scheme, SBL oscillations are described by effective 2-fiavor oscillation prob- 
abilities, involving effective mixing angles for each oscillation channel. The expressions 
for the effective angles 9ee, Gfi/j,, O^ie governing the u e disappearance, u ^ disappearance, 
and — >• V'e appearance probabilities are given in Eqs. (3.2), (4.2), (5.2), respectively. 
From those definitions it is obvious that the three relevant oscillation amplitudes are not 
independent, since they depend only on two independent fundamental parameters, namely 
\Ue4\ and |f7^4|. Neglecting terms of order |?7q,4|^ (a = e,/i) one finds 

sin^ 2(9^e ~ 4 sin^ 2eee sin^ 26*^^ . (6.1) 

Hence, the appearance amplitude relevant for the LSND/MiniBooNE signals is quadrati- 
cally suppressed by the disappearance amplitudes, which both are constrained to be small. 
This leads to the well-known tension between appearance signals and disappearance data 
in the 3+1 scheme, see e.g. [29, 30] for early references. 
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92.4/68 
82.4/68 


7.9 
19.7 
7.8 


616/621 
609/621 
611/621 


10.1 
6.1 
9.0 



Table 7. Global minima, GOF values, and parameter goodness-of-fit (PG) test [125] for the 
consistency of appearance versus disappearance experiments in the 3+1, 3+2, and 1+3+1 schemes. 
The corresponding parameter values at the global best fit points are given in Tab. 8. The last four 
columns give the contributions of appearance and disappearance data to Xpc -'^l- (^•2)- 



This tension is illustrated for the latest global data in the left panel of Fig. 8, where we 
show the allowed region for all appearance experiments (the same as the combined region 
from Fig. 7) , compared to the limit from disappearance experiments in the plane of sin^ 20^e 
and Amii. The preferred values of Am^^ for disappearance data come from the reactor 
and Gallium anomalies. The regions for disappearance data, however, are not closed in this 
projection in the parameter space and include sin^ 2^^e = 4|?7e4f^^4p = 0, which always can 
be achieved by letting — )• because of the non-observation of any positive signal in SBL 
'z^'^ disappearance. The upper bound on sin^ 20^e from disappearance emerges essentially 
as the product of the upper bounds on \Uei\ and \Ufj,4\ from and disappearance 
according to Eq. (6.1). We observe from the plot the clear tension between those data sets, 
with only marginal overlap regions at above 99% CL around Am|^ 0.9 eV^ and at 3a 
around Aml^ w 6 eV^. 

The tension between disappearance and appearance experiments can be quantified by 
using the so-called parameter goodness of fit (PG) test [35, 125]. It is based on the 
definition 

XPG ~ Xmin,glob Xmin,app Xmin,dis AXapp ~l~ AXdis j 

(6.2) 

AXx = Xx.glob - Xmm,x with X = app, dis, 

where Xmin,giob is the minimum of the global data combined, Xmin,app ^nd Xmin.dis 
are the minima of appearance and disappearance data taken separately, and xi giob 
Xx evaluated at the best fit point of the global data. Xpc should be evaluated with the 
number of dof corresponding to the number of parameters in common between appearance 
and disappearance data (2 in the case of 3+1). From the numbers given in Tab. 7 we 
observe that the global 3+1 fit leads to Xmin/dof = 712/680 with a p- value 19%, whereas 
the PG test indicates that appearance and disappearance data are consistent with each 
other only with a p-value of about lO""^. The strong tension in the fit is not reflected in 
the global x^ minimum, since there is a large number of data points not sensitive to the 
tension, which leads to the "dilution" of the GOF value in the global fit, see [125] for a 
discussion. In contrast, the PG test is designed to test the consistency of different parts of 
the global data. 

The conflict between the hints for eV^-scale oscillations and null-result data is also 
illustrated in the right panel of Fig. 8. In red we show the parameter regions indicated by 
the combined hints for oscillations including SBL reactor. Gallium, LSND, and MiniBooNE 
appearance data. Those regions are compared to the constraint emerging from all other 
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Table 8. Parameter values at the global best fit points for the 3+1, 3+2, and 1+3+1 mass schemes. 
7^e is the complex phase relevant for SBL appearance experiments as defined in Eq. (2.2). 




Figure 9. Allowed regions in the plane of |Am4;^| and |Am§]^| in 3+2 (upper-left part) and 1+3+1 
(lower-right part) mass schemes. We minimize over all mixing angles and phases. We show the 
regions for appearance data (light blue) and disappearance data (light green) at 95% CL (2 dof), 
and global data (dark and hght red) at 95% and 99% CL (2 dof). 



data. We find no overlap region at 99% CL. Hence, an explanation of all anomalies within 
the 3+1 scheme is in strong tension with constraints from various null-result experiments. 

6.2 3+2 and 1+3+1 global analyses 

Now we move to the global analysis within a two-sterile neutrino scenario in order to 
investigate whether the additional freedom allows to mitigate the tension in the fit. We 
give s-iicl PG values for the 3+2 and 1+3+1 schemes in Tab. 7 and the corresponding 
values of the parameters in Tab. 8. We observe from the PG values that the tension between 
appearance and disappearance data remains severe, especially for the 3+2 case, with a PG 
value below 10~^, even less than for 3+1. For 1+3+1 consistency at the 2 per mille level 
can be achieved. 

Let us first discuss the 3+2 fit. We find a modest improvement of the total the 
global fit compared to 3+1 by 

X3+i,giob ~ X3+2,giob = 10-''' ■ (6-3) 
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for fixed values of Ar 



Figure 10. Allowed regions for 3+2 in the plane of |C/e4?7^4| vs. 
and Am|]^ at 90% and 99% CL (2 dof). We minimize over all undisplayed mixing parameters. We 
show the regions for appearance data (blue), disappearance data (green), and the global data (red). 



Evaluated for 4 additional parameters relevant for SBL data in 3+2 compared to 3+1 this 
corresponds to 96.9% CL. 

The origin of the very low parameter goodness of fit can be understood by looking at 
the contributions of appearance and disappearance data to Xpc- Tab. 7 shows that the 
of appearance data at the global best fit point, Xappgiob' changes only by about 3 units 
between 3+1 and 3+2. However, if appearance data is fitted alone, an improvement of 
15.2 units in is obtained when going from 3+1 to 3+2, see Eq. (5.3). The fact that 
appearance data by themselves are fitted much better in 3+2 than in 3+1 leads to the large 
value of = 25.8, with a contribution of 19.7 from appearance data. In other words: the 
fit to appearance data at the global 3+2 best fit point (Xappgiob ~ 92.4/68, p- value 2.6%) 
is much worse than at the appearance-only 3+2 best fit point (Xminapp/^of = 72.7/63, 
p- value 19%). This interpretation is also supported by Fig. 6, showing an equally bad fit 
to MiniBooNE neutrino data at the 3+1 and 3+2 global best fit points (black solid and 
red solid histograms, respectively). 

We further investigate the origin of the tension in the 3+2 fit in Figs. 9 and 10. In 
Fig. 9 we show the allowed regions in the multi-dimensional parameter space projected 
onto the plane of the two mass-squared differences for appearance and disappearance data 
separately, as well as the combined region. The 3+2 global best fit point happens close to 
an overlap region of appearance and disappearance data at 95% CL in that plot. However, 
an overlap in the projection does not imply that the multi-dimensional regions overlap. In 
the left panel of Fig. 10 we fix the mass-squared differences to values close to the global 
3+2 best fit point and show allowed regions in the plane of \UeAU^i\ and \Ue^Un^\. These 
are the 5-neutrino analogs to the 4-neutrino SBL amplitude sin20^e- Similar as in the 3+1 
case we observe a tension between appearance and disappearance data, with no overlap 
at 99% CL. This explains the small PG probability at the 3+2 best fit point. The right 
panel of Fig. 10 corresponds to the local minimum of the combined fit visible in Fig. 9 
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Figure 11. Same as Fig. 10 but for tlie 1+3+1 scheme. 

around ^m\^ = 0.9 eV^, Arngj^ = 6 eV^. In this case no tension is visible in the mixing 
parameters shown in Fig. 10, however, from Fig. 9 we see that those values for the mass- 
squared differences are actually not preferred by appearance data, which again leads to a 
degraded GOF. We conclude that the tension between appearance and disappearance data 
cannot be resolved in the 3+2 scheme. 

For the 1+3+1 ordering of 5-neutrino mass states a somewhat better fit can be ob- 
tained. We find 

X3+l,glob ~ Xl+S+l.glob = ^'^■^ 1 (6-4) 

corresponding to disfavoring 3+1 at the 99.9% CL (4 dof) compared to 1+3+1. We observe 
from Tab. 7 that at the 1+3+1 global best fit point a much better fit to appearance data is 
obtained than at the 3+2 best fit point (Xapp giob ~ ^^•'^ compared to 92.4). As visible from 
the blue solid histogram in Fig. 6 the lack of an event excess in the MiniBooNE neutrino 
spectrum around 0.6 GeV is reasonably well reproduced at the 1+3+1 global best fit point, 
although the low energy excess is still under-predicted. The Xpc ^^"^ appearance versus 
disappearance for 1+3+1 is even slightly less than for 3+1 (16.8 versus 18.0). Because of 
the additional parameters relevant for the evaluation of Xpg p-value 0.2% is obtained 
for 1+3+1, about one order of magnitude better than in 3+1. 

The projection of the allowed regions on the plane of the mass-squared differences is 
shown in the lower-right part of Fig. 9. Note that the disappearance regions are to good 
accuracy symmetric for 3+2 and 1+3+1. This can be understood from Eq. (2.3), where the 
difference between 3+2 and 1+3+1 appears only in the last term, which is suppressed by 
the 4th power of small matrix elements, compared to the leading terms at 2nd order. We 
observe in Fig. 9 that appearance and disappearance regions for 1+3+1 both overlap with 
the combined best fit point. In Fig. 11 we show again a section through the parameter 
space at fixed values for the mass-squared differences close to the global best fit point. 
Although the tension between appearance and disappearance is still visible (no overlap 
of the 90% CL regions) the disagreement is clearly less severe than in the 3+2 situation 




-26- 



shown in the left panel of Fig. 10, and in Fig. 11 we find significant overlap at 99% CL, in 
agreement with the somewhat improved PG p-value. 

7 Summary and discussion 

We have investigated in detail the status of hints for eV^-scale neutrino oscillations, namely 
the indications for v e disappearance due to the reactor and Gallium anomalies, and the 
indications for i^^^ — t- z^e appearance from LSND and MiniBooNE. Those hints have been 
analysed in the context of the global data on neutrino oscillations, including short and long- 
baseline accelerator and reactor experiments, as well as atmospheric and solar neutrinos. 
Our main findings can be summarized as follows. 

1. For all fits a global x^^Jdoi « 1 is obtained in our analysis, involving 689 data points 
in total, see table 7. 

2. However, a joint fit of all anomalies suffers from tension between appearance and 
disappearance data, mainly due to the strong constraints from i/^ disappearance 
data. 

3. The tension in the fit is driven by the LSND and MiniBooNE appearance hints, 
since oscillations in the v ^ ^ v channel inevitably predict also a signal in ^ 
disappearance, which is not observed at the relevant LjE scale. 

4. In contrast, the reactor and Gallium anomalies are not in direct conflict with other 
data, since v e and v ^ disappearance at the eV^ scale are controlled by independent 
parameters. 

5. In a 3+1 scheme the compatibility of appearance and disappearance data is at the 
level of 10~^. The individual allowed regions have marginal overlap at about 99% CL. 

6. We do not find a very significant improvement of the fit in a 3+2 scheme compared 
to 3+1. Based on the relative minima, 3+1 is disfavored with respect to 3+2 at 
96.9% CL. The compatibility of appearance and disappearance data in 3+2 is even 
worse than in 3+1, because the fit of appearance data-only is significantly better in 
3+2 than in 3+1, however, the appearance fit at the global best fit point is only 
marginally improved. 

7. We find an improvement of the global fit in the 1+3+1 spectrum compared to 3+1, 
at the 99.9% CL. The compatibility of appearance and disappearance data is still low 
in 1+3+1, at the level of 0.2%. 

Hence, in all cases we find significant tension in the fit, with the marginal exception of the 
1+3+1 scheme. At our 1+3+1 best fit point the minimal value for the sum of all neutrino 
masses would be S 3y^|Am|j^ + y^lAmlj^l + lAm^j^l 3.2 eV, where we took the values 
given in Tab. 8 and assumed that the mass-squared difference with the smaller absolute 
value is negative, using the symmetry 4 5 and 7^/3 — )• — 7^^ of SBL data, see Eqs. (2.1) 
and (2.3). It remains an interesting question whether such a large value of E is consistent 
with cosmology, see e.g. [47, 52, 53, 58, 59]. 
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Let us briefly compare our results to two other recent global sterile neutrino fits, from 
Refs. [34] and [47]. We are in good agreement with the results of [34]. For instance, 
in Tab. 2 of [34] Xpc values for the consistency of appearance and disappearance data 
are given, 17.8 for 3+1 and 23.9 for 3+2, which compare well with our numbers from 
Tab. 7, 18.0 and 25.8, respectively. There is some disagreement with the results of [47]. 
The corresponding Xpc values reported in their Tab. I are 6.6 and 11.12, which lead to 
significantly better compatibility of appearance and disappearance data. Comparing Fig. 1 
of [47] with our Fig. 8 (left) we observe that our disappearance limits are somewhat stronger 
and our appearance region is at somewhat larger mixing angles, both effects increasing the 
tension. Our appearance region is in good agreement with Fig. 6 (left) of [34]. There are 
some differences between our disappearance region and Fig. 6 (right) of [34], mainly at 
high Amii. 

Irrespective of the hints for disappearance and 'i^'^ — )• i^'e appearance, we have 
derived constraints on the mixing of eV-scale states with the r- neutrino flavor. Those are 
dominated by data involving information from neutral-current interactions, which are solar 
neutrino data (NC matter effect and SNO NC data), MINOS long-baseline NC data, and 
atmospheric neutrino data (NC matter effect). The global limit is dominated by the latter. 

In conclusion, establishing sterile neutrinos at the eV-scale would be a major discovery 
of physics beyond the Standard Model. At present a consistent interpretation of all data 
indicating the possible presence of eV-scale neutrino mass states remains difficult. The 
global fit suffers from tension between different data sets. An unambiguous solution to this 
problem is urgently needed. We are looking forward to future data on oscillations at the 
eV^ scale [38], as well as new input from cosmology. 
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A Complex phases in sterile neutrino oscillations 



In this appendix we discuss in some detail the phases for neutrino oscillations involving s 
extra sterile neutrino states. For definiteness, we will focus on s = 2; the special case of 
■s = 1 can be easily obtained by dropping all terms containing a redundant "5" index. 
Let us order the flavor eigenstates as (fg, f/^, i^r, J^sn 1^52) introduce the following 
parametrization for the n x n mixing matrix, with n = 3 + s: 

U = ^35 ^34^25^24^^23^15^44 143^12 (A.l) 

where Vij represents a complex rotation by an angle 9ij and a phase (fij in the ij plane. Note 
that rotations involving only sterile states (i.e., Vw with both £,£' > 4) are unphysical, 
and therefore we have omitted them from Eq. (A.l). Removing those unphysical angles, 
U contains n(n — l)/2 — s{s — l)/2 = 3(s + 1) physical angles. 

In Eq. (A.l) we have chosen a priori all rotations to be complex. We present now a 
method which allows to remove unphysical phases from the mixing matrix in a consistent 
way. First, we note that a complex rotation can be written as 

Vi,=Dk0^jDl, k = i or k=j, (A.2) 

where Oij is a real rotation matrix, Dk is a diagonal matrix with {Dkjjj = e*'^ for j = k 
and {Dkjjj = 1 for j 7^ k. Depending on whether k = i 01 k = j, the phase in is 
either i^pij. Second, we note that phase matrices at the very left or right of the matrix 
U drop out of oscillation probabilities and are therefore unphysical. Hence, we have to 
represent all matrices Vij in Eq. (A.l) using Eq. (A.2), and then try to commute as many 
phase matrices to the left and the right. The matrix Dj. commutes with a matrix Oij if 
k ^ i and k ^ j. Furthermore, ii k = i or k = j we can commute with a complex 
matrix Vij by re-defining the phase fij-. e.g., Vij{6ij,ifij)Di = DiVij{6ij,ip[j). However, we 
cannot commute Dk with a real matrix Oij \i k = i ox k = j. 

This leads to the following rule for removing phases. Let us start by removing one 
phase, let's take for instance ipi2-, obtaining a real V12 — )• Oi2- Then we can no longer use 
the matrices Di and D2 to remove phases, since we cannot commute them with O12 to the 
very left or right of U . But, we can use for instance to remove one of the remaining 
phases (^^3, and so forth. Hence, we can remove in total n — 1 phases. Starting with all 
3(s + 1) physical angles complex, we obtain that there are 3(s + 1) — (n — 1) = 2s + 1 
physical phases, i.e., 1 phase for no sterile neutrinos, 3 phases for the 3+1 spectrum, and 
5 phases for the 3+2 spectrum. Those remaining phases cannot be associated arbitrarily 
to the Vij but only in a way which is consistent with the above prescription to remove 
phases. In particular, it is not possible to make simultaneously three rotation matrices ij, 
ik, kj real. One possible choice is the one given in Eq. (2.6). Using this recipe to remove 
phases it is also straightforward to obtain the physical phases in case of the SBL or LBL 
approximations according to Tab. 2. 

^''in this work we focus on neutrino oscillations. In cases where lepton-number violating processes are 
relevant, such as neutrino-less double beta-decay, more phases will lead to physical consequences and our 
phase counting does not apply. In particular, in such a case the phases on the right of the mixing matrix 
U (these are the so-called Majorana phases) cannot be absorbed. 
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In the SBL approximation for a 3+2 scheme, only two physical phases remain. In 
the parametrization invariant notation from Eqs. (2.1) and (2.2), they are given by 7^e 
and 7^7-- Since the only SBL appearance experiments we consider are studying the 'z^'^ — >• 
V'e oscillation channels only the phase is relevant for our analysis. In the specific 
parametrization from Table 2, the physical phases have been chosen as (^25 and (^35. Since 
(/?35 does not appear in the parametrization independent representation of according to 
Eq. (2.2) we can remove it from our SBL analysis without loss of generality. 

In the LBL limit, more phases are phenomenologically relevant. In particular, Eq. (2.4) 
shows that the oscillation probabilities in the 3+2 case are sensitive to the parametrization 
independent phases 

arg(/a/343 + Ial353) , arg{Ial354:) = 7a/3 , (A.3) 

with la/Sij defined in Eq. (2.2). The experiments for which the LBL approximation is 
relevant are ICARUS and MINOS. ICARUS searches for i/^ — t- transitions, whereas 
the NC data in MINOS are sensitive to the combination Yla=e /i r ^fp-s-i'a • Therefore, for 
our analyses the two appearance channels (a/3) = (/le) and (/xr) are relevant, leading, 
according to Eq. (A.3), to four independent phases, in agreement with Tab. 2.^^ The 
particular parametrization from the table implies that for the i/^ — )• channel only the 
phases (^13 and (^25 are relevant, whereas the i/^ — )• Vr channel is also sensitive to (/J35 and 

^34- 

From the way we have chosen the complex rotations in Tab. 2 the correct phases in the 
3+1 case are automatically obtained by dropping all rotations including the index "5" in 
the 3+2 mixing matrix. We recover the well-known result that in the SBL approximation 
in a 3+1 scenario no complex phase appears. In the LBL approximation two phases remain, 
corresponding to the combinations aig{U*^UeiU^3U*^) and arg{U*^UriU^3U*^) , which can 
be parametrized by using the phases (^534 and (pi3, where for the — >• i^e channel only ipis 
is relevant. 

Let us comment also on the role of phases in solar and atmospheric neutrinos. As 
shown in appendix C solar neutrinos do depend on one effective complex phase. This 
is included in our analysis in full generality however the numerical impact of this phase 
dependence is small. It has been shown in [33] (appendix C) that the impact of complex 
phases on atmospheric neutrinos is very small and we neglect their effect in the current 
analysis. 

B Systematic uncertainties in the reactor analysis 

The correlation of errors between SBL reactors are quite important in order to obtain 
the significance of the reactor anomaly. Here we describe our error prescription for the 
SBLR analysis. From the errors quoted in the original publications we extracted the 
following components. First we removed the uncertainty on the neutrino flux prediction, 
since we include this uncertainty in a correlated way for all reactor experiments based on 

^^In deriving Eq. (2.4) we have assumed that Am^i, Amii, Am§4 are infinite. Note that this assumption 
does not reduce the number of physicai phases further, since also the general procedure used in Tab. 2 
(assuming only Amii = 0) leads to the same number of physical phases as Eq. (2.4). 
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the prescription given in [21] (see below). The remaining error is divided into uncorrelated 
errors (including statistical as well as experimental contributions) as well as correlated 
errors between some SBLR measurements. The total uncorrelated error is shown in the 
last column of Tab. 3. Below we give details on our assumptions on correlations. 

The total error on the measured cross section per fission in Bugey4 is 1.38% [66]. It 
receives contributions which are reactor/site specific (1.09%) as well as detector specific 
(0.84%). Rovno91 [67] used the same detector as Bugey4. The errors on the experimen- 
tal cross section comes from the reactor and geometry (2.1%) and the latter from the 
detector (1.8%). So the first one should be uncorrelated whereas the second one should 
be correlated with the corresponding one from Bugey4. Hence we have o'Bug°y4 = 1.09%, 

^cor n ^uncor o 1 0/ ^cor -i oO/ 

"Bugey4/Rovno91 " U.O-i/o, O^RovnoQl " ^--L/O, O^Rovno91/Bugey4 " ^--O/O. 

The BugeyS measurement consists of 3 detectors at the distances 15 m, 40 m, 95 m. 
In Tab. 9 of [68] systematic errors of 5% (absolute) and 2% (relative) are quoted. The 
uncorrelated errors given in our Tab. 3 are obtained by adding the statistical error (Tab. 10 
of [68]) to the 2% relative systematic error. For the correlated error we remove the relative 
systematic error as well as 2.4% for the fiux prediction and obtain cr^^^^y^ = 3.9%, which 
we take fully correlated between the 3 rate measurements. In cases when we include the 
spectral data from Bugey3 we use 2% (3.9%) as uncorrelated (correlated) normalization 
errors for the three spectra. Details of our spectral analysis of Bugey3 can be found in [79] . 

In Goesgen the same detector was used at three diff'erent distances. In Tab. V of [69] 
the individual and correlated errors are given. The values for the uncorrelated errors 
used in our analysis (see Tab. 3) are obtained by adding the statistical and uncorrelated 
systematic errors in quadrature and expressed in percentage of the ratio. Then [69] quotes 
a correlated error of 6%, which includes 3% from the neutrino spectrum, 2% from the 
cross section, 3.8% from efficiency, 2% from reactor power, and a few more < 1%. We 
remove the 3% neutrino spectrum, as well as the 2% from cross section (this seems way 
too large). This gives CTQ^gggg^^ = 4.8%. Part of this error is supposed to be correlated 
with ILL, since they used a "nearly identical" detector. Removing the reactor power of 2% 
we get fQj^ggggj^yjLL = 4.36%. In the ILL paper [70] errors of 3.66% statistical and 11.5% 
systematical are quoted. The contributions to the systematic error are given as 6.5% on the 
"intensity of the anti-neutrino energy spectrum", 8% detection efficiency, 1.2% neutron life 
time and some other smaller contributions. In the lack of detailed information we proceed 
as follows. We remove 3% for the fiux uncertainties (the same as in Goesgen) and take 
8% (the detection efficiency) to be correlated with Goesgen. This gives cr^-^^°^ = 8.52% 
and (^iij^/Qoesgen ~ '^^ere the uncorrelated error includes also the statistical one. We 
have checked that other "reasonable" assumptions on the ILL / Goesgen correlation do not 
change our results significantly. 

From Krasnoyarsk [71, 72] there are three data points based on a single detector, which 
records events from 2 "identical" reactors. In [71] from 1987, results at distances of 32.8 m 
and 92.3 m are reported. The statistical errors are 3.55% and 19.8%, respectively, and 
the systematical error are 4.84% and 4.76%, respectively, which include detector effects 
(~ 3%), reactor power (~ 3%) and the effective distance (~ 1%). We take systematical 
errors fully correlated between those two data points. Then there is a measurement from 
1994 [72] at 57 m. The errors include detector uncertainty (3.4%), reactor power (2.5%), 
and statistics (0.95%). We assume the detector error to be correlated with the 1987 data 
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points but include the reactor power in the uncorrelated error. 

For SRP [73] measurements at the distances of 18 m and 24 m are reported from the 
same detector, which has been moved between the two positions. The obtained ratios 
of data over expectation at the two distances are 98.7% ± 0.6%(stat.) ± 3.7%(syst.) and 
105.5% ± 1.0%(stat.) lb 3.7%(syst.). The uncorrelated systematic error is derived from the 
ratio of the two spectra, 1.61 ± 0.02(stat) ± 0.03(syst), with an expectation of 1.73 [73]. 
Hence 1.86% = 0.03/1.61 is an uncorrelated systematic error. Then we remove the 2.5% 
from the neutrino spectrum from the systematical error and obtain Cg^pj; = 1.95%, Cgppl = 
2.11%, and Cgf[p = 2.0%. With this assumption on the uncorrelated errors the two data 
points are consistent at about 2. 4a. 

Rovno88 [74] reports 5 measurements with two different detectors: II, 21, IS, 2S, 3S, 
where the "I" experiments use an integral neutron detector, whereas the "S" experiments 
use a scintillation detector measuring the positron spectrum. In Tab. Ill of [74] for each 
measurement two systematical errors are given, 2.2% for "the uncertainty in the measured 
reactor power and the geometric uncertainty" , and a second uncertainty due to "errors in 
the detector characteristics and fluctuations". From Tab. II one finds that statistical errors 
are negligible. In the absence of detailed information we assume the 2.2% uncertainty fully 
correlated among all experiments. From the second error we assume that half of it is 
uncorrelated and the other half is correlated among detectors of the same type. We have 
checked that our results do not depend significantly on those assumptions. 

Finally let us comment on the uncertainty on the neutrino flux predictions. As men- 
tioned above this uncertainty has been removed from the SBLR experimental errors since 
they are treated in a correlated way for all reactor experiments. For the uncertainties of 
the fluxes from ^'^^U, ^^^Pu, ^^^Pu we use the information from tables provided in [21]. 
The uncertainty is provided as uncorrelated error in each bin of neutrino energy as well 
as fully correlated (between energy bins as well as the three isotopes) errors. For the un- 
correlated errors we proceed as follows. We perform a fit of a polynomial of 2nd order 
to the numbers given in [21]. Then those coefficients are used as pulls in the analysis 
constrained by the covariance matrix obtained from the polynomial fits. This allows us to 
take into account the fact that the bin-to-bin uncorrelated errors of the neutrino spectrum 
will lead to correlated effects in the observed positron spectra. Since the uncorrelated flux 
errors are sub-leading compared to the correlated ones the parametrization with a 2nd 
order polynomial is sufficiently accurate. To include the correlated errors we follow [21]: 
the various contributions to this error in each neutrino energy bin are symmetrized and 
added in quadrature. Then we obtain an energy dependent fully correlated error for the 
spectra from ^^^U, ^^^Pu, ^^^Pu which is included as one common pull parameter in the 
global reactor x^- For the neutrinos from ^38^ ^j^^ from [20] and include a 

global normalization error on the 238]j induced events of 8.15% [26]. 

C Solar neutrino analysis 

In the analysis of solar neutrino experiments we include the total rates from the radio 
chemical experiments Chlorine [126], GALLEX/GNO [85] and SAGE [127]. Regarding 
real-time experiments, we include the electron scattering energy-zenith angle spectrum 
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data from all the Super-Kamiokande phases I-IV [128-131] and the data from the three 
phases of SNO [132-134], including the results on the low energy threshold analysis of 
the combined SNO phases I-III [135]. We also include the main set of the 740.7 days of 
Borexino data [136] as well as their high-energy spectrum from 246 live days [137]. In total 
the solar neutrino data used in our analysis consists of 261 data points. 

Let us now focus on the probabilities relevant for the analysis of solar neutrino ex- 
periments. We will assume that only the first two mass eigenstates are dynamical, while 
the others are taken to be infinite. Since physical quantities have to be independent of the 
parameterization of the mixing matrix, we will use the freedom in choosing a parameteriza- 
tion that makes analytical expressions particularly simple. We start from the Hamiltonian 
in the flavor basis: 

H = UAU^ + V , (C.l) 

where A = diag(0, Am^^, Amj^, . . .)/2E and V = \/2 diag(2iVe, 0, 0, Nn, . . . )/2. It 
is convenient to write U = UU12, where U12 is a complex rotation by an angle 612 and a 
phase 612 which we will define later. Then we can write: 

H = UHW with H = UuAul^ + WVU . (C.2) 

In order to further simplify the analysis, let us now assume that all the mass-squared differ- 
ences involving the "heavy" states Uh with h > 3 can be considered as infinite: Am^^ — )• 00 
and Am^^, — t- 00 for any I = 1,2 and h, h' > 3. In leading order, the matrix H takes the 
effective block-diagonal form: 

'H'-'l 



where H^'^'^ is the 2x2 sub-matrix of H corresponding to the first and second neutrino 
states, and A*^'^^ = (Ma.g{Am\i, Am|]^, . . . )/2E is a diagonal (s + 1) x (s + 1) matrix (the 
matter terms in this block are negligible in the limit of very large Am|^, ) . Consequently, 
the evolution matrix is: 



S ^ _ .^(.)^ with 5(2) = Evol [^(2)] and S = USW . (C.4) 



We are interested only in the elements Sae- It is convenient to define 9i2 in such a way 
that Ue2 = 0. Taking into account the block-diagonal form of S, we obtain: 

5„e = f/el (UalS?j^ + Ua2Si^A + ^ ^/^/.L^.^e"*^'^'^^ . (C.5) 



The expressions for the probabilities, Pae = 1 5*0.6 1 , are straightforward: 



h>3 

2 



Pae — l^ell 



Ual\^\S{^^\^ + \Ua2\^\Sil^\^ + 2Re {UalUa2S[^^ Sif*^)] + ^ \Uah\^\Ueh\' 

h>3 

Uel\' [{\Ua2\' - |^.l|2)|5(?|2 + 2Re(?7,i^,*25fi^5(f )] + ^ \Ua^\W. 



i=all 

(C.6) 



'^^We define O12 in such a way that the corresponding rotation is the rightmost rotation in U, but we do 
not assume anything about aU other rotations. If the standard parameterization is used for them, then 612 
has the usual interpretation. 
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Here we have used the fact that the terms containing a factor e^*^'^'^^ oscillate very fast, 
and therefore vanish once the finite energy resolution of the detector is taken into account. 
For solar neutrino experiments we only need Pee and Pae = Pee + Pfie + Pre- It is therefore 
convenient to define 612 in such a way that Yl(7=st ^0-1^^2 ^ ^^^^ number. Let us also 
define: 

Ve = \Uelf, ^D=Y1 (l^<^2|' - Palf) , CiV = 2 ^ U^iU*^ . (C.7) 

lT=St <T=st 

Using unitarity relations, |C/mp = 1 and ^„ UaiU*2 = 0, we obtain: 

Pee = Ce - r,lP^} , Pae = Ca - r,e {^dP^'J + CnP^^) (C.8) 

where 



i=all i=all iT=st 



issions 

H^'^\ which has the form: 



(C.9) 



(2) (2) 

In the above expressions Pose and P;„4. are effective terms derived from the Hamiltonian 



IT=St 



a2 \Ua2\' 



(C.IO) 



with the vacuum term including the phase Si2- 

(2) ^ Am|i /- cos 2^12 sin 29^2 e'^^^ \ 

AE V sin 2^12 6-^-512 cos 2^12 ; ^ ^ ^ 

These equations arc valid for any number of extra sterile states. Hence, in the most general 
case solar neutrinos depend on six real mixing parameters, 

^12, IJe, ^D, Cn, Ce, Ca , (C12) 



one complex phase (5i2), and Am2i. Let us now further define Vs ^ y S,d + S,n ^^'^ discuss 
some special limits: 



?7e = 1: this implies Uei = 1, hence Uei = for i > 3. By unitarity |f7ai| = 1) 

hence Uai = for a 7^ e. Therefore, Ce = Ca = 1 and = 0; the probabilities 

(2) (2) 
reduce to the well-known expressions Pee = 1 — Pose and Pae = 1 — ??sPosc with 

Vs = ^d; 

• ris = 0: in this case = = 0, sterile neutrinos do not participate to 1-2 oscillations 
and Pae is a constant. If we further assume that no sterile neutrino exists then 
Ce = ?7e + (1 — ?7e)^ and Ca = 1, so that Pae = 1 and Pee reduces to the well-known 
three-neutrino formula; 
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• s = 1: for one sterile neutrino r]s reduces to rjg = \Usi\'^ + |t^s2p- 

To perform the analysis we map the parameters which are used for the analysis into 
the effective parameters for solar neutrinos. In the case of 3+1 the number of real mix- 
ing parameters is actually the same as the number of mixing angles in the most general 
parametrization of the 4x4 mixing matrix: the six parameters in Eq. (C.12) are a function 
of the six angles 9i2, 6*13, ^14, O23, O24, The dependence of solar neutrinos on 6*13,^14 
is shown in Fig. 3 and the one on 024, ^34 follows from Fig. 5. The dependence on 623 is 
important for the NC matter effect and SNO NC data. For s > 1 the number of effective 
mixing parameters of solar data is less than the number of angles in the general mixing 
matrix. The phase ^12 is a complicated function of complex phases and angles. We have 
verified numerically for the 3+1 case that if all angles are non-zero (and ^24, ^34 relatively 
large) the from solar data varies by about 1 to 2 units as a function of the phase. Once 
all relevant constraints on the mixing angles are imposed the effect of the phase on solar 
data is negligible. See also [138] for a discussion of complex phases in solar neutrinos in 
the context of sterile neutrinos. 



D Atmospheric neutrino analysis 

The analysis of atmospheric data follows closely the one presented in Refs. [139, 140] and 
includes the Super-Kamiokande results from phases I, II and III [141] (80 data points in 
total). Technical details on our fit can be found in the appendix of [142]. 

In order to derive suitable expressions for the relevant probabilities, we can follow the 
approach presented in App. C for solar data. The Hamiltonian in the flavor basis is given 
by Eq. (C.l), which in the mass basis becomes: 

H = UHU'< with H = A + U'fVU . (D.l) 

As before, we can simplify the analysis by assuming that all the mass-squared differences 
involving the "heavy" states Vh with h > A can be considered as infinite: Arnf^- — )• 00 and 
Am^^, —7- c« for any i = 1,2,3 and h,h' > 4. In leading order, the matrix H takes the 
effective block-diagonal form: 

where H^^^ is the 3x3 sub-matrix of H corresponding to the first, second and third neutrino 
states, and A'^'^-' = diag(Am|]^, Am^i, . . . ) /2£'^ is a diagonal sxs matrix (the matter terms 
in this block are negligible in the limit of very large Am^^, ) . We are interested only in the 
probabilities P^/j with a,/3 G {e,//,r}. Taking into account the block-diagonal form of H, 
we obtain: 



Paf} = \Aap\^ + K^p , ^ = ?7(3)^(3)c/(3)t , K^f} = ^\U^h\^\Uph\^ (D.3) 

/i>4 



where 5*^3) = Evol [^(3)j ^nd C/^^) the 3x3 sub- matrix corresponding to the first three 
lines and columns of U, and as such it is not a unitary matrix. 
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Let us now focus on the three-neutrino system described by H^^^ = A^^^ + V^^\ with 
= [C/t\/f/j(3)_ 'pj^g matter term V^^^ contains both the "standard" contribution from 
z/e charged-current interactions, and a "non-standard" part induced by the absence of 
neutral-current interactions for sterile neutrinos. In practical terms, this is just the same 
as the problem of neutrino propagation in the presence of non-standard neutrino-matter 
interactions (NSI) described in Ref. [140]. Following the approach discussed there, we make 
a number of simplifying assumptions: 

• we assume that the neutron-to-electron density ratio, Rne, is constant all over the 
Earth. We set Rne = 1-051 as inferred from the PREM model [143]; 

• we set Am^i = 0, thus forcing the vacuum term A^^^ to have two degenerate eigen- 
values; 

• we impose that the matter term V^^^ also has two degenerate eigenvalues. 

The first two assumptions are very well known hand have been discussed in detail in the 
literature. For example, in Sec. 5.2 of Ref. [142] it was noted that the different chemical 
composition of the Earth mantle and core has very little impact on NSI results. The 
last approximation is adopted here for purely practical reasons, since (together with the 
other two) it allows to greatly simplify the calculation of the neutrino evolution [140, 144]. 
Unfortunately, in the present context the matter term V^^'^ cannot be fixed a priori, but it 
arises as an effective quantity determined by the angles and phases of the mixing matrix 
U. Finding the points in the general parameter space for which V^^'^ has two degenerate 
eigenvalues (the only points for which our numerical analysis is technically feasible) is not 
an easy task. To solve this problem, we have considered here two alternative cases, both 
based on physically motivated scenarios: 

(a) decouple the electron flavor from the evolution and include the NC matter effect; 

(6) allow the electron flavor to participate in oscillations but neglect the NC matter 
effect. 

Option (a) requires to set Uei = for z > 3 (in addition to Am-gi = 0), whereas option 
(6) is equivalent to setting Rne = ("hydrogen-Earth" model). Both approximations have 
been previously discussed in App. C of Ref. [33]. Although we know that none of these 
options corresponds to Nature, we can make a sensible choice of when it is safe to use each 
of them. It turns out that for constraining the mixing of the i'^ with eV-scale neutrinos 
the NC matter effect plays no role at all (see discussion in App. C2 of [33]), whereas 
the participation of the electron flavor may have some impact. Therefore, whenever we 
are mainly interested in constraining |C/^/i| {h > 4), as in the case of the global analysis 
combined with SBL data, we adopt assumption (6). This is important since non-zero \Uei\ 
leads to slightly relaxed constraints on |^7^4|, although the effect is small once external 
constraints on \Uei\ are taken into account. With respect to our former analysis presented 

^^The reason is that the main information comes from effects on the overah normahzation of /i-hke 
events, and in order to fix this normahzation the e-hke sample plays a role in constraining uncertainties in 
the neutrino fluxes. 
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in Ref. [33], the explicit NSI formalism adopted here is more general since it allows to fully 
include |C/ei I "'^slated effects without further approximations. 

On the other hand, when exploring the sensitivity to the fraction of sterile neutrinos 
participating in atmospheric and long-baseline neutrino oscillations, the contribution of 
neutral-current neutrino-matter interactions is essential. Indeed, under approximation (6) 
no limit on \Uri\ {i > 4) would be obtained. Therefore, when exploring constraints on \Uri\ 
we adopt assumption (a) above. This is relevant for Figs. 4 (right) and 5. 

E Technical details on the simulation of SBL and LBL experiments 

Here we provide technical details on the simulation of some of the experiments included in 
our fit. All the simulations described in this appendix make use of the GLoBES software 
package [145, 146]. 

E.l KARMEN/LSND Ue disappearance analysis 

Both LSND and KARMEN have measured the reaction i^e + ^^C ^ + ^^N, where the 
^^N decays back to ^^C + e'^ + i^e with a lifetime of 15.9 ms. By detecting the electron from 
the first reaction, which has a Q-value of 17.33 MeV, one can infer the neutrino energy. 
Here we describe how we use these data to constrain disappearance [94, 95]. 

For KARMEN, we use the information from the thesis [94], which uses more exposure 
than the original publication [92]. The number of expected events can be calculated by 
multiplying the ^^C cross section (Fukugita et al. [93]: 9.2ibl.l x 10"^'^ cm^), the number of 
target nuclei (2.54x 10^°), the absolute neutrino flux (5.23x lO^i), the efficiency (27.2%, fiat 
in energy), and the inverse effective scaled area (l/[47r(17.72m)^]). 846 neutrino candidates 
are observed, with an expected background of 13.9 it 0.7, which are mainly accidentals and 
cosmic induced. The systematic errors are dominated by a 6.7% uncertainty in the absolute 
neutrino flux and 3% in the Monte Carlo efficiency. The total systematic error is 7.5% plus 
a 12% cross section error. We take the latter to be correlated between the KARMEN 
and LSND i/e-carbon analyses. In fig. 3.2 (upper panel) of [94], the data are shown as the 
energy distribution of the prompt e~ spectrum in 26 bins where the visible prompt electron 
energy is within the range 10 MeV < Eg < 36 MeV. Modulo energy reconstruction the 
neutrino and electron energies would be related by the Q-value as = Eg + Q- For the 
simulation we assume 30 MeV < E,^ < 56 MeV. To properly fit the data, we assume 
a"e = 1h%l\jE (MeV) for the energy resolution. With the 26 data points we obtain a 
two-neutrino best fit with Xmin/^of = 30/24. 

For LSND [91] we compute the expected number of events by multiplying the ^^C cross 
section (9.2 x 10^^^ cm^ [93]), the number of target nuclei (3.34 x 10'^''), the neutrino fiux 
at the detector (10.58 x 10^^ cm-^), and the efficiency (23.2%, flat in energy). 733 neutrino 
candidates are observed, with a negligible expected background. The systematic error is 
dominated by a 7% uncertainty in the neutrino flux and a 6% uncertainty in the effective 
fiducial volume. The total systematic error, not including the theoretical cross section 
error, is 9.9%. The 12% cross section error is correlated between the LSND and KARMEN 
fe-^^C analyses. In fig. 6 of Ref. [91], the data are shown as the energy distribution of 
the prompt e~ spectrum, where the visible prompt electron energy is within the range 
18 MeV < Ee < ^2 MeV and divided into 12 bins of width 2 MeV. In terms of the neutrino 
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energy E^, this energy range corresponds to 35.3 MeV < Ey < 59.3 MeV. In our analysis, 
we combine the 12 energy bins into only 6 bins. To properly fit the data, we assume 
(Te = 2.7 MeV for the energy resolution. With the 6 data points we obtain a two-neutrino 
best fit with Xmin/dof = 3.81/4. 

For the combined KARMEN+LSND i/g-carbon fit, we have 32 bins and we find a 
two- neutrino best fit point with Xmin/*^°f ~ 34.17/30. 

E.2 E776 

The pion beam experiment E776 at Brookhaven [40] employed a 230 ton calorimeter de- 
tector located at approximately 1 km from the end of the 50 m long pion decay pipe. The 

V e energy of order GeV was measured with an energy resolution of 20% / E [GeV] . E776 
used *z/'^ disappearance data to obtain the overall normalization of the neutrino fiux. In 
our fit, we do not explicitly include 'z/'^ data, but instead use the normalization as an 
input. The main backgrounds in E776 came from intrinsic z/'e contamination in the beam 
and from vr'^'s produced in neutral current interactions and misidentified as electrons. The 
systematic errors were 11% for the intrinsic background and 27% (39%) for the vr'^ back- 
ground in neutrino (anti-neutrino) mode. In 1986, E776 collected 1.43 x 10^^ (1.55 x 10^^) 
protons on target, and a total of 136 (56) Ve (i^e) candidate events were observed with an 
expected background of 131 (62) events for neutrino (anti-neutrino) mode. E776 present 
the observed and predicted electron energy spectra using 14 equidistant energy bins per 
polarity, covering the energy range from GeV to 7 GeV. In our fit we omit the first bin 
and combine the second and third ones because modeling the detection efficiency at these 
low energies is very difficult. Hence, we have a total of 24 data points. We checked that 
we are able to reproduce well the exclusion curve shown in fig. 4 of ref. [40] (if we also use 
a two-flavor oscillation model), and we obtain Xmin/dof = 31.08/22 at the best fit point. 
In the combined analysis with other experiments we take into account oscillations for the 

V e background. 

E.3 ICARUS 

The ICARUS experiment [147] is a neutrino beam experiment at Gran Sasso. The CNGS 
facility at CERN shoots 400 GeV protons at a graphite or beryllium target, producing 
a hadronic shower which is focused by a magnetic horn system. The resulting neutrino 
beam is mainly composed of f^, having only a 2% z/^ contamination and a < 1% intrinsic 
Ve component. The neutrino spectrum ranges approximately from — 50 GeV, with a 
wide peak at 10 — 30 GeV. After traveling 732 km, they are detected in the ICARUS 
T600 detector, a 760 ton liquid argon time projection chamber. Between 2010 and 2012, 
the ICARUS detector observed 839 neutrino events with energy below 30 GeV, to be 
compared with the expectation of 627 and 3 Vr charged current events, as well as 204 
neutral current events. While a Monte Carlo simulation of the experiment predicts 3.7 Uf. 
background events, only two were identified [41] . For our ICARUS simulation, we took the 
Vfj, spectrum from Ref. [148] and folded it with the — )• z/g oscillation probability. To assess 
the limit on v^. appearance, we define the likelihood — 21n(L) = 2{P — D) + 2Dlog{D/P), 
where D (P) is the total observed (expected) number of fe events. Although there is a 
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~ 7% systematic error in the selection efficiency, we have checked that its impact on the 
experiment sensitivity is negligible. 

E.4 MiniBooNE 

To analyze MiniBooNE data on — oscillations, we use the data presented in [16], 
corresponding to 6.46 x 10^^ protons on target in neutrino mode and 11.27 x lO^'^ protons 
on target in anti-neutrino mode, the same exposure as used in [17]. We follow the analysis 
strategy outlined in the supporting on-line documentation [16]. For each set of oscillation 
parameters we consider, we compute the expected neutrino and anti-neutrino energy spec- 
tra by weighting the unoscillated Monte Carlo events from MiniBooNE's data release [16] 
with the oscillation probabilities. For each simulated event, we take into account the in- 
dividual neutrino energy and the distance between the neutrino production and detection 
vertices. We allow both signal and background neutrinos to oscillate including the "wrong- 
sign contamination" of the signal. In computing the log-likelihood for any given parameter 
point, we take into account statistical and systematic uncertainties by using the covariance 
matrix published by MiniBooNE [16]. The likelihood is given by 

-2ln{L) = (D - Pf S-\n - P) , (E.l) 

where D (P) is a column vector containing the observed (predicted) number of events 
including background, and S is the covariance matrix (see Ref. [16] for details). Note that 
this covariance matrix includes correlations between the and z^'^ event samples as well 
as between neutrino and anti- neutrino data (if analysed together). There are 11 bins for 
e-like data and 8 bins for /i-like data, for neutrino and anti-neutrino mode, each. 

The I' ^-data are used to determine the normalization of the beam flux. In principle, 
we should therefore include also oscillations in the i^'^ disappearance sample. The im- 
pact of i^'^ oscillation on the fit to the z^'g data would then be taken into account by the 
covariance matrix. This is problematic because it would prevent us from combining, with- 
out double-counting, our results with the independently obtained ones from MiniBooNE's 
dedicated i^^^ disappearance searches (see below). Using the u ^ data from [16] directly 
for a disappearance analysis is also not possible because the corresponding prediction 
from [16] has been obtained from a Monte Carlo simulation in which some parameters have 
been tuned to the data assuming no disappearance. We have therefore decided to fol- 
low the MiniBooNE collaboration and not to include oscillations of in the appearance 
analysis. We have verified numerically that this has a negligible impact on the fit results 
once external limits on |C/^4|, |C^/j5| are taken into account. 

A further subtlety arises for the kaon-induced backgrounds. In MiniBooNE, these 
are predicted using information from the SciBooNE detector [149], which operated in the 
same beam, but at a shorter baseline. To account for this, we rescale the kaon-induced 
backgrounds in MiniBooNE by the ratio of the oscillation probabilities in MiniBooNE 
and in SciBooNE. Since kaon-induced backgrounds make only a small contribution to 
MiniBooNE's total error budget, the effect of this rescaling on the fit results is, however, 
very small. 

Due to the correlation with the muon-like events it is not straight forward to assign a 
number of degrees of freedom to the appearance search without double counting the muon 



-39- 



data, which are used also in the separate disappearance analysis. We have adopted the 
following prescription. Eq. (E.l) can be written as 

= -2 ln(L) = d^Meede + 2dlMe^,d^ + d^M^^d^, 

= {de + 6YMee{de + 5) + C 

where de {d^) are the e (/x) components of the vector (D — P), M = S'"^, and M^p denotes 
the corresponding sub-blocks of the matrix M. In Eq. (E.2) we have defined 

5 = M-^^Mef,d^ and C = dj^ {M^^ - M^,M-^Me^,) d^ = d^iS^^y^d^ , (E.3) 

where {S^^)~^ is the inverse of the jj.^ sub-block of SM Hence, we have block-diagonalized 
the covariance matrix. The shift 6 corresponds to the impact of the //-like data on the 
normalization of the e-like flux. The two terms in Eq. (E.2) should be statistically in- 
dependent and approximately distributed. For the MiniBooNE appearance analysis 
we therefore use Xmb app = ~ C, and assign 22 dof to it (for combined neutrino/anti- 
neutrino data). The last equality in Eq. (E.3) shows explicitly that C does not depend on 
oscillation parameters, since we neglect the effect of oscillations on d^. With this method 
we obtain GOF values which are in reasonable agreement with the numbers obtained by 
the collaboration: our results for Xmin/'^o^ (GOF) for neutrino, anti- neutrino, combined 
data are 14.2/9 (11%), 6.5/9 (69%), 32.9/20 (3.5%), respectively, compared to the numbers 
obtained in [16] 13.2/6.8 (6.1%), 4.8/6.9 (67.5%), 24.7/15.6 (6.7%). Note that in [16] the 
number of dof and GOF have been determined by explicit Monte Carlo study and also a 
different energy range has been used to obtain those numbers. 

For our MiniBooNE disappearance analysis, we use the neutrino mode data from [44]. 
As in appearance mode, we compute the expected event spectra for each parameter point 
by using MiniBooNE's Monte Carlo events. Since backgrounds are very small for this 
analysis, we do not need to take them into account. For each set of oscillation parameters, 
we choose the overall normalization of the spectrum in such a way that the total predicted 
number of events matches the number of observed events in MiniBooNE, i.e., we fit only 
the event spectrum, not the normalization. The log-likelihood is obtained in analogy to 
Eq. (E.l) and thus takes into account systematic uncertainties and correlations between 
different energy bins. 

In the analysis of i'^ disappearance data, we follow the combined MiniBooNE/SciBooNE 
analysis from [45]. As for the appearance and disappearance analyses, we use public 
Monte Carlo data to compute the predicted event spectra. We take into account oscilla- 
tions of both the signal and the background, and we compute the log-likelihood again in 
analogy to Eq. (E.l). 

E.5 MINOS 

Our analysis of MINOS neutral current (NC) and charged current (CC) interaction is based 
on 7.2 X 10^'' protons on target of NC data presented by the collaboration in Ref. [150] 
(see also [42, 43]) and 7.25 x 10^° protons on target of CC data published in [151]. All 

'^^Note that (5*^^)^^ is difFerent from M^^, the latter being the fifi sub-block of the inverse of S. 
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data was recorded in neutrino mode, i.e., the beam consists mostly of and only small 
contaminations of i/^, I'e and Ue- 

We have implemented the properties of the NuMI beam and the MINOS detector 
within the GLoBES framework [145, 146], using results from the full MINOS Monte Carlo 
simulations as input wherever possible. In particular, we use tabulated Monte Carlo 
events [152] to construct the detector response functions R^^{E^'^^'^, E^'^'^) and R^^ {E^'^'^'^ , E^'^^) 
for neutral current events in the near detector (ND) and the far detector (FD), respectively, 
i?^^ and i?^^ describe the probability for a neutrino with true energy E^'^^'^ to yield an 
event with reconstructed energy E'^'^^°. We include the CC and CC backgrounds to 
the NC event sample (beam intrinsic as well as oscillation induced), as well as the small 
NC background to the CC event sample. The number of charged current interactions is 
predicted using the simulated NuMI flux [153], the cross sections calculated in [154-156], 
and a Gaussian energy resolution function with width ag^/E*''^^ = 0.1/ Y^^;*''"^/GeV for 
the CC event sample and = 0.16 + 0.07/ v^£;*™VGeV for the CC background 

in the NC event sample. The parameters of the energy resolution function, as well as the 
efficiencies, have been tuned in order to optimally reproduce the unoscillated event rates 
predicted by the MINOS Monte Carlo. We emulate the baseline uncertainty due to the 
non-vanishing length of the MINOS decay pipe by smearing the oscillation probabilities 
with an additional Gaussian of width cj^ = (2.0 GeV — 3.0 x E^'^^'^)'^ and setting the eff'ective 
distance between near detector and neutrino source to 700 m. This value as well as the 
parameters of the smearing function have been obtained numerically from the requirement 
that 2-flavor oscillation probabilities at the near detector computed with a more accurate 
treatment of the decay pipe are well reproduced for various mass squared differences of 
order eV. 

We compute neutrino oscillations in MINOS numerically using a full 4- or 5-flavor 
code that includes all relevant mixing parameters as well as CC and NC matter effects. 
In the fit, we predict the expected number of events at the far detector for a given set of 
oscillation parameters by multiplying the observed number of near detector events in each 
energy bin with the simulated ratio of far and near detector events in that bin. Each event 
sample is divided into 20 energy bins with a width of 1 GeV each, covering the energy range 
from to 20 GeV. As systematic uncertainties, we include in the analysis of the NC (CC) 
sample a 4% (10%) overall normalization uncertainty on the far-to-near ratio, separate 15% 
(20%) uncertainties in the background normalization at the far and near detectors, a 3% 
(5%) error on the energy calibration for signal events and a 1% (5%) error on the energy 
calibration for background events. For the NC analysis, the systematic uncertainties are 
based on the information given in [42], for the CC analysis they have been tuned in order 
to reproduce the collaboration's fit with reasonable accuracy, while still remaining very 
conservative. 
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